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ABSTRACT 


This  report  examines  the  problem  of  adaptively  tracking  a 
maneuvering  submarine  in  two  dimensional  space  utilizing  passive  time 
delay  and  Doppler  frequency  measurements  of  unknown  or  randomly 
varying  center  frequencies.  The  target  is  free  to  maneuver  in 
velocity  and  depth  with  tracking  being  done  in  the  vertical  plane. 

It  is  pointed  out  how  to  incorporate  bearing  measurements  into  the 
present  polar  model  to  achieve  a  three  dimensional  target  tracking 
capability. 
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Chapter  1 


1.1  Introduction 

During  the  past  several  years  much  effort  has  been  spent  in  the 
development  of  sophisticated  digital  filtering  algorithms  for  tracking 
maneuvering  targets.  A  common  method  has  been  to  model  the  target 
dynamics  in  a  rectangular  coordinate  system  which  results  in  a  linear 
set  of  state  equations,  but  forces  the  measurements  to  be  nonlinear 
functions  of  the  state  variables.  With  this  model  an  extended  Kalman 
filtering  algorithm  is  frequently  used  both  to  provide  current  state 
variable  estimates  and,  by  a  one-step  prediction  process,  to  linearize 
the  next  measurement  vector.  This  method  works  moderately  well  until 
the  target  makes  an  abrupt  change  in  its  trajectory  in  response  to 
pilot  or  missile-guidance  program  commands.  In  this  situation  the 
velocity  and  position  estimates  can,  and  often  do,  diverge  from  the 
true  unknown  values.  The  inherent  problems  of  this  approach  can  lead 
to  large  bias  errors  and  sometimes  complete  filter  divergence. 

Earlier  work  on  the  maneuvering  target  tracking  problem  includes 
Jazwinski’s  limited  memory  filtering  [1],  in  which  the  filter  gains 
are  prevented  from  decaying  to  zero.  Another  technique,  described 
by  Thorp  [2],  involves  switching  between  two  Kalman  filters  in  res¬ 
ponse  to  a  detected  maneuver.  A  third  approach,  due  to  Singer  [3], 
models  the  target  trajectory  as  a  response  of  the  target  model  to 
a  time-correlated  random  acceleration.  With  this  method  additional 
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state  variables  are  used  to  generate  the  correlated  forcing  func¬ 
tions  which,  in  turn,  increase  the  dimension  of  the  Kalman  filtering 
algorithm.  In  this  manner  the  technique  provides  the  filter  with 
statistical  information  concerning  target  maneuvers  based  on  an  as¬ 
sumed  range  of  possible  aaelerations.  Singer's  method  was  subsequen¬ 
tly  extended  by  many  others. 

Parallel  to  the  effort  was  the  method  of  modeling  major  changes 
in  target  trajectories  by  a  semi-Markov  process.  An  application  of 
this  approach  to  tracking  maneuvering  targets  in  two-dimensions  by 
Moose  [4]  was  successfully  extended  by  Gholson  and  Moose  [5]  to 
three-dimensional  tracking . 

The  general  approach  which  uses  the  "adaptive  semi-Markov  man¬ 
euver  model"  of  [4]  and  [5]  implies  a  discretization  of  possible 
vehicle  accelerations  or  velocities.  The  estimation  algorithm  then 
views  the  maneuvering  vehicle  as  if  it  is  responding  to  commands  which 
are  modeled  by  a  semi-Markov  process,  i.e.,  a  random  process  with  a 
finite  number  of  "states"  (commands)  which  are  selected  according  to 
the  transition  probabilities  of  a  Markov  process.  A  semi-Markov  pro¬ 
cess  differs  from  a  Markov  process  in  that  the  duration  of  time  in 
one  state  prior  to  switching  to  another  state  is  itself  a  random  var- 
ble  [6].  Incorporating  the  semi-Markov  concept  into  a  Bayesian  esti¬ 
mator  of  [4]  and  [5].  This  estimation  algorithm  provides  a  substan¬ 
tial  improvement  in  filter  stability,  which  means  that  large  bias 
errors  are  prevented  from  being  built  up  due  to  unmodeled  target  accel¬ 
erations.  An  important  aspect  of  this  adaptive  estimation  algorithm 
is  its  elimination  of  a  "growing  memory"  which  is  prevalent  in  many 
adaptive  filters. 
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1.2  Target  Modeling 

With  the  brief  history  of  the  maneuvering  target  tracking  pro¬ 
blem  presented  in  the  previous  section,  we  see  a  general  progression 
in  the  sophistication  of  tracking  filter  design  stemming  primarily 
from  the  method  in  which  the  unknown  target  accelerations  are  modeled. 
This  trend  is  graphically  outlined  in  Figure  1.1 

Initially,  target  maneuvers  were  modeled  as  the  response  to 
uncorrelated,  zero-mean  variations  about  a  nonaccelerating  target, 
shown  in  Figure  1(A).  As  a  result,  the  estimation  algorithm  could 
follow  only  those  maneuvers  which  were  comparable  with  the  input  noise 
level.  Furthermore,  the  filtering  results  during  nonmaneuvering  sit¬ 
uations  were  seriously  degraded  due  to  the  uncorrelated  input  noise. 

As  shown  in  Figure  1(B),  Singer  [3]  attempted  to  model  large-scale 
maneuvers  by  assuming  a  time-correlated  input  process  and  incorpora¬ 
ting  the  statistics  into  the  subsequent  filter  design.  In  Figure  1(C) 
large-scale  target  maneuvers  were  modeled  as  a  stochastic  process  whose 
mean-value  switched  randomly  from  among  a  finite  set  of  predetermined 
values.  The  adaptive  estimation  algorithm  mentioned  in  the  previous 
section  could  then  be  used  to  track  the  maneuvering  target.  This 
method  was  seriously  restricted,  however,  by  the  requirement  of  a 
large  number  of  preselected  mean  values  in  order  to  ensure  convergence 
of  the  estimation  process.  In  this  report  we  show  that  by  combining 
the  concepts  illustrated  in  Figure  1(B)  and  Figure  1(C)  the  number  of 
mean  values  required  to  prevent  filter  divergence  is  greatly  reduced. 
This  combination  is  illustrated  in  Figure  1(D).  The  primary  benefit 
of  this  approach  is  the  large  saving  in  computational  effort.  An 
additional  benefit,  at  least  from  a  subjective  viewpoint,  is  that  the 
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Figure  1.1.  Historical  development  of  maneuvering  target  model  inputs. 
(A)  Zero-mean  white  Gaussian  plant  disturbance.  (B)  Correlated,  zero- 
mean  plant  disturbance.  (C)  White  Gaussian  noise  with  randomly  switching 
mean.  (D)  Correlated  Gaussian  noise  with  randomly  switching  mean. 


Figure  1.2.  Target  motion  model. 


time-correlated, randomly  switching.mean-forcing  function  more  ade¬ 
quately  models  real-world  target  maneuvers. 

The  basic  target  modeling  ideas  are  shown  in  Figure  1.2.  The  tar¬ 
get  trajectory  is  generated  by  the  random  selection  of  an  input  time- 
correlated  Gaussian  process  whose  mean  value  is  applied  to  the  tar¬ 
get  plant  dynamics  for  a  random  duration  of  time.  This  input  distur¬ 
bance  process  lasts  until  a  new  input  Uj  is  randomly  chosen  from  among 
a  finite  set  of  n  possible  inputs.  With  this  model  as  a  background 
and  using  an  appropriate  choice  of  state  variable  equations  to  repre¬ 
sent  target  dynamics,  either  submarine  or  aircraft,  it  is  possible 
to  develop  an  "optimal"  (in  the  minimum  mean-square  error  sense) 
tracking  filter  that  adaptively  learns,  then  quickly  adjusts  itself 
for  each  major  alteration  of  target  trajectory. 


1.3  Incorporation  of  Singer  Process  into  the  Target  Dynamics 
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is  a  white  Gaussian  random  process  acting  in  the  x  direction 


A  similar  set  of  equations  exists  for  the  y  and  z  directions. 
Def ining 
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the  following  continuous  time  state  variable  model  is  obtained  for 
Equation  (1.3.1) 
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Discretizing  (1.3.2)  in  time  yields 
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where 


A  *  (1  -  e"aT)/a 


B  -  [1  +  (ae  aT  -  ae  aT)  /  (a  -  a) ]  /  (aa) 


(aT  -  1  +  e_aT)/a2 


D  ■  [T  +  (aA  -  oJ)/(a  -  a)]/(aa) 


G  -  (J  -  A)/a  -  a) 


Generation  of  a  Correlated  Gaussian  Random  Process  w'  by  Filtering 
the  White  Gaussian  Random  Process  w 


Chapter  2 


UNDERWATER  PASSIVE  TRACKING.  INTRODUCTION  AND  LINEARIZED  POLAR  MODEL 
FOR  THE  MANEUVERING  OBSERVER- SOURCE  SCENARIO 

2.1  Introduction 

The  tracking  of  air  targets  such  as  described  in  [7]  is  charac¬ 
terized  by  high  quality  radar  data  providing  good  system  observability 
and  high  signal-to-noise  ratio.  In  underwater  passive  submarine  trac¬ 
king  on  the  other  hand,  the  situation  is  characterized  by  meager  data 
which,  when  it  is  available,  is  almost  always  a  nonlinear  function  of 
the  system  state  variables . 

The  particular  problem  which  is  addressed  in  this  chapter  and 
in  the  remainder  of  this  report  is  how  to  track  the  target  submarine 
(Source)  from  the  moving  Observer  submarine  using  sonar  signals  eman¬ 
ating  from  the  Source .  Because  the  Observer  does  not  use  an  active 
sonar  but  rather  listens  to  sound  waves  emanating  from  the  Source, 
this  method  of  tracking  has  acquired  the  name  passive  tracking. 

It  is  pointed  out  in  [7]  in  conjunction  with  air  targets  that  in 
order  to  successfully  track  a  maneuvering  target  in  three  dimensions 
the  dynamics  must  be  realistically  modeled  in  the  chosen  coordinate 
system  and  secondly  that  an  adaptive  feature  must  be  designed  to  pre¬ 
vent  filter  divergence  as  the  target  executes  evasive  maneuvers.  These 
two  requirements  are  of  equal  importance  in  the  underwater  passive  trac¬ 
king  domain. 
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The  availability  of  air  target  radar  data  in  spherical  coordi¬ 
nates  suggests  that  coordinate  system  as  being  the  "natural"  one  in 
which  to  design  a  tracking  filter.  In  the  underwater  case,  the  selec¬ 
tion  of  a  "natural"  coordinate  system  based  on  the  available  data  is  a 
much  more  obscure  process .  This  is  true  because  almost  all  the  avail¬ 
able  data  are  nonlinear  functions  of  the  state  variables  regardless  of 
any  coordinate  system  choesn.  It  is  a  major  goal  of  this  report  to 
address  this  problem  in  detail.  It  will  be  shown  that  for  a  certain 
subset  of  the  observables,  a  "natural"  coordinate  system  does  indeed 
exist.  This  discovery  will  then  be  exploited  to  design  a  "pseudo- 
linear"  filter  to  process  this  subset  of  the  observables  using  only 
linear  filtering  techniques. 

2.2  Stare  of  the  Art 

T 

The  problem  of  passively  tracking  underwater  targets  has  re¬ 
ceived  very  little  attention  in  the  literature.  What  little  exposure 
has  occurred  has  tended  to  concenetrate  on  the  bearings-only  approach. 
In  this  method  the  Observer  monitors  his  bearing  to  the  Source,  over 
a  period  of  time.  Usually  the  Observer  must  execute  a  known  maneu¬ 
ver  to  remove  the  ambiguity  that  sometimes  exists  with  this  type  of 
tracking  algorithm. 

In  a  very  interesting  paper,  Hassab  [8]  passively  tracks  the 
maneuvering  Source  in  three  dimensions  by  exploiting  the  combination  of 
bearing  measurements  and  sonar  time  delays  between  the  direct,  bottom 
reflected  and  surface  reflected  sound  waves  emitted  by  a  Source.  Since 
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these  time  delays  and  bearing  measurements  are  nonlinear  functions  of 
the  state  variables  in  the  rectangular  coordinate  system  used.  The 
Extended  Kalman  Filter  is  employed  to  provide  optimal  estimates  of 
the  Source  range,  bearing,  velocity  and  depth  with  respect  to  the  moving 
'Observer.  The  state  vector  consists  of  three  rectangular  position  com¬ 
ponents  (R  ,  R  ,  R  )  and  three  rectangular  velocity  components  (V  ,  V  , 
x  y  2  x  y 

V  ).  A  constant  velocity  model  is  used  for  the  Source  dynamics  with 
an  additive  white  Gaussian  noise  perturbation.  Because  of  the  choice 
of  a  rectangular  coordinate  system,  the  processing  of  the  bearing  mea¬ 
surements  becomes  coupled  to  the  processing  of  the  time  delay  measure¬ 
ments  resulting  in  a  3  x  6  linearized  measurement  matrix  having  the 
following  generic  form: 


y 


0  0 


H(k+i) 


x 


0 


0 


(2.2.1) 


where  3,  t^,  are  the  bearing  and  two  sonar  time  delays,  respectively. 

In  [9]  a  relative  rectangular  coordinate  system  is  again  chosen 
to  model  the  Observer-Source  scenario.  The  measurements  used  here  are 
bearing  and  direct/surface  reflected  sonar  time  delays.  Since  these 
make  the  system  unobservable,  the  Source  is  assumed  to  be  at  the  sane 
depth  as  the  Observer.  Again,  the  measurements  are  nonlinear 
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functions  of  the  state  variables  with  coupling  occurring  between  the 
processing  of  bearing  and  time  delay  measurements. 


One  of  the  factors  which  makes  rectangular  coordinates  appear  so 
attractive  for  passive  tracking  is  that  compensation  for  Observer  motion 
in  this  coordinate  system  is  relatively  simple,  involving  only  a  linear 
subtraction.  For  example,  if 

V  »x  component  of  Source  velocity  with  respect  to 
s  the  ocean  floor 

V  *  x  component  of  Observer  velocity  with  respect  to 
o  the  ocean  floor 

R  ■  relative  distance  between  the  x  coordinates  of  the 
Source  and  Observer 

then,  assuming  a  constant  velocity  model  for  the  Source  [8] 


n* 

p  — 

r  — 

““ 

R  (k+l) 

1  T 

R  (k) 

T 

X 

3 

X 

- 

V  (k+l) 

0  1 

V  (k) 

0 

Lxs  J 

-  — 1 

L  xs  - 

V  (k) 
x 


(2.2.2) 


Thus  rectangular  coordinates  offer  a  very  convenient  method  of 
compensation  for  Observer  motion. 

However,  rectangular  coordinates  are  a  very  poor  coordinate  sys¬ 
tem  in  which  to  filter  when  one  considers  the  types  of  measurement  data 
other  than  sonar  time  delays  available  for  passive  tracking.  For  ex¬ 
ample,  Source  bearing  with  respect  to  the  Observer  is  one  of  the  most 
consistently  available  measurements.  Another  commonly  available  mea¬ 
surement  data  consists  of  Source  line  frequency  spectra  which  are 
emitted  from  various  mechanical  elements  aboard  the  Source. 

In  filtering  in  relative  rectangular  coordinates  (x,  y,  z) 
these  types  of  measurement  data  are  inherently  nonlinear  functions  of 
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the  Source  state  variables.  For  example,  bearing  involves  the  arc 
tangent  function  of  Source  x  and  y  relative  coordinates.  The  line  fre¬ 
quency  spectra,  when  used  in  Doppler  tracking  with  an  unknown  center  fre¬ 
quency  fQ,  involve  the  produot  of  state  variables  and  again  a  nonlinear 
estimation  problem  must  be  solved. 

In  an  attempt  to  find  a  more  suitable  coordinate  system  for 
passive  tracking,  Tenney  et  al.  [10]  propose  a  relative  coordinate  sys¬ 
tem  based  on  the  closest  point  of  approach  (CPA)  for  use  with  Doppler 
tracking.  This  approach  requires  that 

1.  the  measurement  matrix  be  relinearized  at 
each  iteration; 

2.  the  Source  trajectory  be  of  the  "crossing" 
variety  to  produce  a  closest  point  of 
approach; 

3.  the  center  frequency  of  the  transmitted 
signal  be  known  and  constant. 

The  presence  of  bearing  observations  suggests  either  polar  or 
spherical  coordinates  as  a  more  suitable  system  in  which  to  filter. 

This  is  true  in  light  of  earlier  work  where  a  linearized  spherical  fil¬ 
ter  was  used  to  process  spherical  measurement  data  using  linear  filtering 
theory  while  permitting  decoupling  of  each  coordinate  direction.  Thus, 
for  example,  the  bearing  measurements  available  to  the  Observer  could 
be  decoupled  from  the  time  delay  measurements  and  processed  independ¬ 
ently  in  either  of  these  two  coordinate  systems.  While  the  sonar  time 
delay  measurements  would  continue  to  be  coupled  together  and  necessitate 
a  nonlinear  filter  such  as  the  Extended  Kalman  filter,  the  decoupling 
of  the  bearing  tracker  from  the  time  delays  would  represent  a  signi¬ 
ficant  reduction  in  the  complexity  of  the  nonlinear  filter.  For 
example,  if  filtering  in  polar  coordinates  is  undertaken,  then  because 
of  the  decoupling  of  the  bearing  channel  the  time  delay  processing 
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would  involve  linearizing  the  measurement  matrix  about  only  two  state 
variables  rather  than  about  three  state  variables  as  in  Equation  (1.2.1). 

To  distinguish  between  the  merits  of  spherical  and  polar  co¬ 
ordinates,  consider  Figure  2.3.1.  Assume  that  both  the  Observer  and 

Source  depths  remain  unchanged;  therefore  z  9  0,  If  the  Source  radial 

so 

•  « 
rate  o  is  non-zero  then  in  spherical  coordinates  the  elevation  rate  e 

is  also  non-zero,  even  if  p  is  some  fixed  value.  In  polar  coordin¬ 
ates,  on  the  other  hand,  and  with  these  same  conditions,  the  Source  ver¬ 
tical  race  is,  of  course,  zero.  The  implications  of  all  this 
are  thac  in  spherical  coordinates  the  generally  non-zero  elevation  rate 
requires  a  sophisticated  modeling  technique  for  the  Source  input 
acting  in  the  elevation  direction.  This  technique  would  involve  several 
overlapping  Gaussian  curves  and  also  an  adaptive  technique  to  provide  a 
weighted  sum.  It  will  be  shown  in  Chapter  3  that  with  polar  coordinates 
on  the  ocher  hand,  the  source  Z  direction  input  can  be  modeled  in  a  par¬ 
ticularly  simple  manner  using  only  a  single  Gaussian  curve  having  a 
mean  value  equal  to  zero. 

It  will  also  be  shown  later  chat  polar  coordinates  possess  ex¬ 
tremely  attractive  features  for  passive  tracking  using  the  Doppler 

effect.  For  these  reasons  it  is  decided  to  undertake  passive  tracking 
of  a  moving  Source  from  a  moving  Observer  in  the  polar  coordinate  sys¬ 
tem. 

Of  course,  while  polar  coordinates  offer  the  advantages  out¬ 
lined  above,  there  is  a  concomitant  penalty  in  that  compensation  for  Ob¬ 
server  motion  in  this  coordinate  system  is  more  difficult  and  complex 
than  for  the  rectangular  case.  To  appreciate  this  point,  consider  the 


state  model  given  by  Equation  (2.2.2).  If  one  imagines  a  rectangular 
coordinate  system  attached  to  and  moving  with  the  Observer,  then  the 
state  model  (2.2.2)  provides  estimates  of  the  target  state  variables 
with  respect  to  this  moving  coordinate  system.  However,  regardless  of 
the  maneuvers  executed  by  the  Source  and  Observer,  this  coordinate  sys¬ 
tem  always  translates  parallel  to  itself. 

Referring  to  Figure  2.3.1,  consider  a  polar  coordinate  system 
attached  to  and  moving  with  the  Observer.  As  the  Observer  and  Source 
maneuver  with  respect  to  each  other  the  direction  of  the  polar  radius 
vector  p  connecting  the  Observer  to  the  Source  changes  its  direction  in 
keeping  with  the  changing  relative  positions.  Therefore  the  polar  co¬ 
ordinate  system  attached  to  the  Observer  rotates  in  addition  to  trans¬ 
lating  and  it  is  the  added  complexity  caused  by  this  rotation  that  ren¬ 
ders  more  difficult  the  task  of  formulating  a  state  variable  model  in 
polar  coordinates  which  accurately  compensates  for  Observer  motion. 

It  is  this  problem  which  is  solved  in  the  next  section. 

2.3  Linearized  Polar  Model  for  the  Dynamics  of  the  Maneuvering 
Observer-Source  Scenario 

Consider  Figure  2.3.1  which  shows  the  geometry  of  the  Observer- 
Source  scenario  in  polar  coordinates.  In  this  figure,  the  plane  defined 
by  the  X-Y  axes  is  parallel  to  the  ocean  bed  and  fixed  with  respect 
thereto.  The  horizontal  distance  separating  the  Source  and  Observer  is 
labeled  p  and  will  be  referred  to  in  the  future  simply  as  the  radius. 
This  polar  radius  is  to  be  distinguished  from  the  spherical  radius 
labeled  r  which  is  the  distance  separating  the  Observer  and  Source  in 
three-dimensional  space.  The  vertical  distance  zgQ  is  simply  the  diff¬ 
erences  in  their  depths.  Figure  2.3.2  is  the  projection  of  Figure  2.3.1 
onto  the  X-Y  plane.  The  parameters  appearing  in  this  figure  are  defined 


in  Table  1.3.1. 
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Figure  2.3.1  Geometry  of  the  Observer-Source  Scenario  in  Polar  Coordinates 


The  linearized  polar  model  is  developed  using  a  modified  ver¬ 
sion  of  the  approximate  modeling  technique  of  Reference  (5)  except  that 


now,  the  origin  of  Che  coordinate  system  is  moving. 


P  CHANNEL  MODEL 


p,  calculation: 
k+1 

Referring  to  Figure  1.3.2 


P  =* 


So 


CCx,  -  x  )2 
s  o 


(ys  -  y0> 


(ys 


,2,1/2 

yj  3 


So 

Sx 

o 


p 


(2.3.1) 
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Figure  2.3.2  Geometry  of  the  Observer-Source  Scenario  Projected  Onto  the  X-Y  Plane 


Table  2.3.1 


Explanation  of  Parameters  Appearing  in  Figure  3-3.2 


p  the  polar  radius 

Vq  the  velocity  vector  of  the  Observer 

3  the  angle  between  the  polar  radius  p 

and  the  (fixed)  X  axis 


SQ  angle  between  the  X  axis  and  the 

velocity  vector  of  the  Observer 

Sgo  the  relative  bearing  of  the  Source 

with  respect  to  the  Observer 

an  x-directed  vector,  used  for 
illustrative  purposes 

V  a  y-directed  vector,  used  for 
illustrative  purposes 

ig  a  unit  vector  in  the  direction  of 

increasing  angle  3 

V  projection  of  the  Observer  velocity 

8  vector  onto  ig 


position  coordinates  of  the  Observer 
and  Source,  respectively,  in  the 
rectangular  coordinate  system 


(x  -  x  )  (y  -  y  ) 

P  m  —  - -  (x  -  x  )  +  -  (y  -  v  ) 

P  s  o  p  ws 

is  expanded  as  follows  keeping  only  Che  linear  terms: 


(2 .3.2) 


k+1 


Pk  +  hT\  (xs  “  xs  )  +  !r"l 
K  *Vk  Sk+1  sk 


T - !  vv  -  V  ) 

3Vk  sk+l  ‘sk 


(2.3.3) 


Assuming  a  linear  drag  model  for  che  Source  as  given  in  (1.3.1),  upon 
substituting  (2.3.1)  into  (2.3.3)  for  che  Source  connected  terms,  and  we  get 
<\-H  -  V  '  '“k  +  3“xk  +  C% 


Sk+1  '  °k  * 


(x  -  X  ) 

- 1  [Ax  +  Bw'  +  Cu  +  Dw  1 

P  'k  s  s  s  s  J  i, 

X  x  x  'k 


(ys  -  y  ) 

+  - - - j  [Ay  +  Bw'  +  Cu  +  Dw  ] 


3  S 

y  y 


(x  -  X  ) 

■  — - —  L  [x  -  X  ] - - 

s  k  °k+l  V 


y  1  k 

(yc  -  yj 


£_!  £yn  -  y„  1 

k  k+l  °k 


(2-3.4) 

In  che  above  equation  and  in  all  subsequent  analysis,  subscripts  s  and 
o  refer  to  the  Source  and  the  Observer,  respectively. 


Now 


(x  -  x  )  *  x  T 

O,  ,  O,  0, 

x+1  k  ec 


and 


(2-3.5) 


(y„ 


y  )  *  v  T 

1  n  >  •  Q 

k 


'k+1  k 

Combining  terms  with  like  coefficients  in  (2.3.4)  and  making  use  of 
(2.3.5), 
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p 


(x  “X 

P,  , .  ■  0.+  A[ — - —  x  + 

lc+1  k  1  p  3 


(x  -  O  (y  -  y  ) 

+  B[— 2 - —  w'  +  — 2 - 2_  w»  ] 

O  a  a  * 


(x  -  X  )  (y  -  y  ) 

+  G  —  u  +—1— — £-u  ], 

P  Sx  °  sy  ]k 


(x3  "  Xo)  (y3  "  V 

+  D[ - - -  V  +  -  W  ]  | 

P  3  0  S  1  , 

x  y  k 


<**  -  V  .  .  -  v  .  , 

_  - - -  x  +  — - -  y  ] 

P  o  p  7oJ 


(2.3.6) 

Consider  che  coefficient  of  B  appearing  in  Equation  (2.3.6).  From  Fig 

ure  2.3,2,  ^X£.  X°^  and  y°^  are  che  direction  cosines  between 

P  P 

che  X  and  Y  axes,  respectively,  and  che  o  direction.  Hence 


(XS  ’  X0)  W’  +  (7s  -  V  W- 


(2.3.7) 


is  the  sum  of  che  projections  of  w'  and  w'  onto  the  p  direction.  Th 

3x  3y 

sum  acting  in  the  o  direction  can  be  replaced  by  a  single  equivalent 
term  denoced  by  .  In  a  similar  manner  the  coefficients  of  C  and  D 

p 

are  called  u  and  w  respectively.  The  coefficient  of  A  in  (2.3.6) 
s.  s 

P  P 

represents  che  projection  of  che  Source  velocity  onto  the  radial  direc¬ 
tion.  This  coefficient  must  be  recast  in  a  different  formulation  in 
order  to  complete  che  state  model  of  (2.3.6).  To  this  end  Equation 


(2.3.2)  is  rewritten  as  follows: 


21 


(x  -  x  )  (y  -  y  ) 

using  — 2 -  i  + - ■ — ■  y  5  V  Cos  3 

p  o  p  oo  so 


(2  .3.8) 

(2.3.9) 


using  2.3.8  and  2.3.9  yields 

ills - Hsi  x  +  — - —  V  =*  p  +  V  CosB  (2.3.10) 

p  s  p  '  s  o  so 

Substituting  Equations  (2-3.9)  and  (2.3.10)  for  the  coefficients  of 

T  and  A,  respectively,  in  (2.3.6)  and  collecting  like  terms  results  in 

the  following  state  variable  model: 

9k+l  "  Pk  +  +  BWs  +  Cus  +  ^s  +  CA"T1(V0Cos8so)  'k  (2,3,U) 

Pk  Pk  Pk 


A  similar  approach  is  taken  to  develop  a  state  model  for  o^+^ 

By  assuming  that  the  Observer  maintains  a  constant  velocity  for  long 
periods  of  time,  and  Chat  is  a  zero  mean  correlated  GAUSSIAN  random 
process  acting  in  the  p  direction  the  following  state  equations  are 
easily  obtained. 


"  p  - 

"1  A  B  “ 

-p 

“  C  (A-T)* 

f  Us  1 

D 

SP 

p 

* 

0  E  F 

o 

A  (E-l) 

V  CosS 

G 

o  so 

k 

w' 

-« 

1 

o 

o 

w' 

0  0 

J 

s_ 

s 

_  * 

— 

k+1 

L.  __ 

_  p„ 

k 

* 

(  2.3.12) 

The  underlying  constraint  for  the  state  model  (2.3.12>  requires  that 
the  Observer  adopt  a  constant  velocity  profile. 
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Z  Channel 


Using  a  discretized  version  of  the  basic  linearized  drag  model 
of  chapter  1,  and  letting  zgQ  be  the  relative  separation  between 
Source  and  Observer  Che  following  state  model  is  presented. 


(2.3.13) 

The  state  model  is  also  subject  to  the  constraint  of  a  constant 
velocity  Observer. 

It  is  interesting  to  observe  that  both  state  models  (2.3.12) 
and  (2.3.13)  are  identical  although  the  types  of  manipulations  involved 
in  the  derivations  are  quite  different.  For  example,  a  Taylor  series 
expansion  is  used  to  derive  (2.3.12)  whereas  no  such  expansion  is 
used  in  deriving  (2.3.13). 

Bearing  Channel 

From  a  similar  development  the  resultant  state  equations  are  presented. 


1 


1 


3  so 

"l 

A 

<B/pkr 

CO 

OJ 

°  a 

(C/pk} 

(A-T) 

“so 

m 

0 

E 

(F/pk) 

3so 

+ 

(A/pk) 

(E-l) 

w’ 

L  3s  J 

k+1 

0 

0 

e~aT 

w' 

s3 

k 

0 

0 

V  Sing 

(B  -  -2 - 2% 

o  p 


Jk 


+ 


"(D/pk) 

(G/pk) 

J 


(2.3.14) 


This  ststs  mod  si  is 


This  completes  the  derivation  of  the  linearized  polar  model  in  relative 
coordinates  for  the  moving  Observer-Source  scenario.  In  the  next  chap¬ 
ter  this  model  will  be  applied  to  passively  track  a  maneuvering  Source 
using  sonar  time  delays . 


2.4  Conclusion 

Before  closing,  a  comment  needs  to  be  made  concerning  the  con¬ 
straints  which  the  preceding  state  variable  models  place  on  the  Observer 
motion.  While  the  state  variable  models  are  only  valid  for  a  non¬ 
maneuvering,  i.e.  constant  velocity  Observer,  this  does  not  mean  that 
the  Observer  cannot  execute  any  maneuvers.  What  is  implied  is  that 
during  the  maneuver  the  state  models  are  inaccurate  and  will  not  yield 
good  estimates.  Upon  completion  of  the  Observer  maneuver,  the  model  is 


23 


Il 


once  again  valid  and  will  yield  good  estimates.  Of  course,  during  a 
maneuver  transients  are  introduced  and  consequently  any  filter  would 
inevitably  yield  poor  estimates  during  this  period. 


! 


Chapter  3 

PASSIVE  TRACKING  USING  SONAR  TIME  DELAYS.  THE  ADAPTIVE. 
EXTENDED  POLAR  KALMAN  FILTER. 


3.1  Introduction 

In  this  chapter  the  linear  polar  model  which  was  developed  In 

Chapter  2  will  be  employed  to  passively  track  a  maneuvering  Source  using 

sonar  time  delays.  Figure  3.1.1,  which  uses  notation  similar  to  that 

in  Hassab  [  3 ] ,  gives  the  two-dimensional  geometry  showing  the  paths 

traversed  by  the  Source-emitted  sonar  signals  in  question.  The  "direct" 

signal  traverses  the  path  labeled  r.  The  "surface  reflected"  signal 

traverses  the  path  labeled  R„  -  R.  and  t.  is  the  difference  in  time  of 

2s  Is  1 

arrival  between  this  signal  and  the  "direct"  signal,  as  measured  at  the 
Observer.  The  "bottom  reflected"  signal  traverses  the  path  labeled 
R2b  ~  Rlb’  an<*  t2  *S  C^e  difference  in  of  arrival  between  this 

signal  and  the  direct  signal,  as  measured  at  the  Observer. 

From  Hassab's  work  [  8]  the  equations,  which  are  repeated  below 
for  convenience,  show  that  and  are  nonlinear  functions  of  the 
system  state  variables  p  and  zgo  as  given  in  Figure  3.1.1. 
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Traversed  by  the  Source-Emitted  Sonar  Signals 
in  Reaching  the  Observer 
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Therefore  the  Extended  Kalman  Filter  must  be  employed  to  provide  optimal 
estimates  in  the  presence  of  these  nonlinear  measurements. 

The  '^adaptive  filter  approach  developed  in  [  5]  was  earlier 
applied  with  much  success  to  tracking  maneuvering  air  targets.  This 
approach  will  now  be  applied  in  the  underwater  environment  to  passively 
track  a  maneuvering  Source  from  a  moving  Observer  using  the  measurements 
(3.1.1)  and  (3.1.2).  Because  of  the  coupling  in  (3.1.1)  and  (3.1.2)  between 

the  p  and  2go  state  variables,  the  p  and  z  channels  cannot  be  filtered 
independently.  Bather  the  two  models  developed  separately  in  Chapter  2 
must  be  combined  in  the  manner  described  in  Section  3.2. 

A  modified  version  of  the  conventional  Extended  Kalman  Filter 
algorithm  must  be  developed  in  order  to  accommodate  nonlinear  filtering 
using  the  adaptive  filter  approach.  This  will  also  be  described  in 
the  next  section. 

3.2  Combined  p/2  Channel  Filtering  of  Sonar  Time  Delays  Using  the 
Augmented  Extended  Polar  Kalman  Filter 

Consider  the  following  combined  p/z  state  vector: 

•  •  T 

[p  p  w'  z  z  w'  ] 
s  s 

P  z 

where  the  subscript  so  on  zgQ  has  been  dropped  for  notational  conve¬ 
nience.  The  linearized  state  variable  model  corresponding  to  this 
state  vector  is  given  in  Equation  (3.2.1).  Defining 

-  $x£i}  +  ru£i}  +  Ywk 

where : 
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Consider  the  deterministic  input  vector  (3.2.1).  The  first 

entry,  namely  u  refers  to  the  mean  value  of  the  ith  Gaussian 

SP 

curve  in  the  series  of  N  partially  overlapping  Gaussian  curves 

now  being  used  to  model  the  Source  p  channel  input.  Likewise  the 

superscript  i  on  ug  ^  in  (3.2.1)  refers  to  the  ith  Gaussian 
z 

curve  which  is  being  used  to  model  the  Source  z  channel  input. 

The  last  entry,  namely  zq,  is  the  (known)  Observer  velocity  in 
the  Z  direction  or  vertical  direction. 

In  general,  the  Source  does  not  execute  many  maneuvers  in 

the  Z  direction.  Whenever  such  a  maneuver  does  take  place  it  lasts  only 

for  a  brief  period.  This  is  true  because  the  Source  can  move  upward 

only  as  far  as  the  ocean  surface  and  downward  to  its  maximum  permissible 

depth,  which  is  usually  a  relatively  small  distance  compared  to  the 

ocean  depth.  Therefore  the  Source  is  quite  constrained  in  its  movement 

along  the  Z  direction  in  Figure  3.1.1.  The  conclusion  to  be  drawn  from 

this  is  that  the  Source  input  ug  is  generally  zero  with  there  being 

z 

only  brief  periods  of  time  when  it  is  non-zero.  Therefore  in  modeling 
the  Source  input  in  the  Z  direction,  only  one  time  correlated  Gaussian 
curve  with  a  mean  value  of  zero  is  needed.  However,  since  the  measure¬ 
ments  are  nonlineup  functions  of  the  state  variables  the  extended 
Kalamon  filter  is  used. 
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3.3  Summary  of  the  Extended  Kalman  Filter 

The  Extended  Kalman  filter  is  usually  the  first  filter  to  be 
used  when  confronted  with  either  system  or  measurement  nonlinearities, 
or  both.  This  filter  is  well  documented  [11]  and  is  summarized 
in  (3.2.2)  -  (3.2.7)  for  a  linear  system  with  nonlinear  measurements. 


P(k+l|k+l)  -  (I-K{X(k+l| k) }H{X(k+l jk) }]P(k+l | k) 

X  [I-K{X(k+l|k)}H{X(k+l|k)}]T 

4-  K{X(k+ljk)}R(k+l)K{X(k+l|k)}T  (3.2.2) 


K(X(k+l|k)}«  P(k+1 1 k)HT{X(k+l J k) } [H{X(k+l jk) }P(k+l | k)H{X(k+l | k) }T 


+  R(k+1) ]  x  (3.2.3) 

T 

P(k+ljk)  -  «P(kjk)*‘  +  'fQ(k)’?  +  (3.2.4) 

X(k+l|k+l)  -  X(k+ljk)  +  K(X(k+l|k)}[z(k+l)-h{X(k+ljk)}]  (3.2.5) 

X(k+1 | k)  «  *X(k|k)  +  Fu(k)  (3.2.6) 

3h  (X) 

HU(k«|k)}i  (-1_1|-(k+l|k)  (3.2.7) 


Several  comments  need  to  be  made  concerning  the  above  algorithm: 

1.  Equation  (3.2.2)  for  the  updated  covariance 
involves  more  terms  than  are  usually  used, 
namely 

P(k+l|k+l)  -  [I-K{X(k+l | k) }H{X(k+l j k) } ]P(k+l jk)  (3.2.8) 

The  additional  terms  appearing  in  (3.2.2)  which 
do  not  appear  in  (3.2.8)  are  necessary  co 
guarantee  that  P(k+l|k+l)  is  always  symmetric. 

The  loss  of  symmetry  that  may  occur  in  using 
(3.2.8)  in  the  Extended  algorithm  arises 
from  the  approximate  nature  of  this  filtering 
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algorithm,  especially  the  approximate 
value  for  the  Kalman  Gain  matrix  which 
the  algorithm  is  only  capable  of  pro¬ 
ducing. 


2.  The  Kalman  Gain  matrix  is  shown  to  be  a 
function  of  the  predicted  estimate 
X(k+l|k)  by  explicitly  including  this 
term  in  its  argument.  This  functional 
relationship  arises  from  the  presence 
in  Equation  (3.2.3)  of  the  linearized 
measurement,  matrix  H{X(k+l|k)}.  This 
linearized  measurement  matrix  is  defined 
by 


3h. (X) 

H{X(k+l| k, }  i  t-5^-]lt(k+1|k) 


3.  The  measurement  residuals  are  formed  by 
evaluating  the  nonlinear  functional 
relationships  in  (3.2.5)  using  the  pre¬ 
dicted  estimate. 


4.  The  additional  term  FQurT  in  (3.2.4)  is 
uaeu  LO  account  for  Cue  mismatch  of 
u(i)  and  the  actual  unknown  u. 


With  the  above  comments  complete,  the  adaptive  Extended  Kalman  filter 
algorithm  is  now  discussed. 

3.4  The  Adaptive  Tracking  Filter 

The  heart  of  the  adaptive  filter  developed  in  this  report  is 
the  forming  of  the  estimate  of  the  target  states,  in  each  channel,  from 
a  weighted  sum  of  estimates  conditioned  on  the  N  possible  discrete  input 
levels. 


Consider  the  state  model  (3.2.1).  This  state  model  views  the 
target  input  acting  in  the  polar  direction  as  being  derived  from  a 
time  correlated  Gaussian  density  having  a  mean  value  up .  Next  consider 
a  series  of  N  such  Gaussian  curves  with  displaced  mean  values 
u  i  *  1,  2,  ...»  N  and  partially  overlapping  "tails"  as  shown  in 

Figure  3.2.1.  If  a  bank  of  N  Kalman  filters  is  formed,  each  filter 
based  on  the  state  equations  of  equation  (3.2.1)  with  the  deterministic 
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input  ^  being  a  different  one  of  these  N  mean  values,  then  a  series 
of  N  estimates  is  obtained,  each  conditioned  on  a  different  Gaussian 
curve.  Next  a  weighted  sum  of  these  estimates  is  obtained  in  a  manner 
to  be  disclosed  below,  and  this  weighted  sum  is  taken  to  be  the 
unconditioned  estimate  of  the  target  states. 

Calculation  of  weighting  coefficients:  The  general  case 

As  previously  discussed,  the  target  input  is  now  modeled  as 
coming  from  a  series  of  N  overlapping  Gaussian  curves  each  of  which 
has  a  predetermined  mean  value.  As  the  target  executes  a  series  of 
evasive  maneuvers  in  the  polar  channel,  for  example,  the  changing  input 
to  produce  these  maneuvers  is  viewed  as  randomly  switching  among  these 
N  curves.  By  applying  the  semi-Markov  statistics  to  this  switching 
process  a  series  of  N  probabilities  W^,  i  =  1,  2,  .  .  . ,  «  is  generated 
[  7 ]  where 

S  Pr  {target  input  is  being  derived  from  the  Gaussin 

(i) 

curve  whose  mean  value  is  u  . } 

P 

These  are  then  used  to  form  the  weighted  estimate. 

We  begin  with  the  well-known  relation  that  the  optimal  estimate 

can  be  written  as  a  weighted  sum  of  the  input-conditioned  estimates 

.  A  (i) 

as  shown  in  figure  3.2.2  [113.  Thus  if  X  (k+1)  represents  the  optimal 
estimate  of  X(k+1)  given  that  the  ith  input  force  u^  is  present 
(the  ith-input  force  being  one  of  the  previously  described  mean  values), 
then  based  on  the  data  sequence 

Z(k+1)  =>  {z(l),  z(2) ,  ....  z(k),  z(k+l)}, 
we  define 

N 

x(k+l)  *  l  xU;(k+l)  W  (k+1)  (3.2.9) 

i-1  1 
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where 


(3.2.10) 


W  (k+1)  -  Pr  (u(k)  -  u(i)|Z(k+l)} 

and 

x(i)(k+l)  -  E{x(k+1) | u(k)  -  u(i),  Z(k+1) } . 

Equation  (3.2.9)  is  a  total  probability  expression  developed  from  the 
basic  relation  that 

x(k+l)  =  E{x(k+1) | Z(k+1) } 

is  the  optimal  mean-squared  estimate.  It  is  well  known  that  the 
optimal  input-conditioned  estimates  are  provided  by  suitably  matched 
Kalman  filters.  In  particular,  for  the  i—  filter 

x(i)(k+l)  -  4>(k)x(i)(k)  +  T(k)u(l)  +  K(k+1)  [z(k+l) 

-  H4(k)x(i)(k)  -  Hr(k)  u(l)]  . 

where 

K(k+1)  =  M(k+1)  HT[HM(k+l)  HT  +  R]"1, 

M(k+1)  »  <t(k)  P(k)  $T(k)  +  'P(k)  Q¥T(k) 

and 

P(k+1)  -  [I  -  K(k+1)H]  M(k+1) . 

The  matrices  T  and  ¥  are  used  to  denote  the  respective  coefficient 
matrices  in  (3.2.1). 

The  following  is  an  outline  of  the  analysis  given  in  [  5  ]  to 
calculate  the  recursive  weighting  coefficients  Wi#  i  =  1,  2,  ...,  N. 
Defining  Z(k+1)  *  (Z(k),  z(k+l)},  apply  Bayes  Theorem  to  (3.2.10)  and 


obtain 


m 


The  denominator,  which  varies  with  time;  is  independent  of  i  and  is 
therefore  common  to  each  W^(k+1)  as  a  normalizing  constant.  The  first 
numerator  factor  is  determined  from  the  semi-Markov  input  process. 
Expanding  this  factor  in  a  total  probability  expression, 

N 

Pr{u(k)  «  u(i)|z(k)}  -  l  Pr{u(k)=u(i) |u(k  -  l)=u(j),  Z(k)}W.(k). 

j-1  2 

And  since  Z(k)  has  no  influence  on  the  Markov  state  transitions. 


Pr{u(k)  -  u(i)|Z(k)}  =  l  6  W.(k) 

j-1  2  2 


(3.2.12) 


where  the  semi-Markov  probability  is 


9  *  Pr{u(k)  *  u^|u(k-l)  =  u^} 


Combining  (3.2.11)  and  (3.2.12) 


W  (k+1)  =  C  p{z(K+i)  |  u(k)  =  ua>,  Z(k)}  £  3  W  (k)  (3.2.13) 

j-1  2  2 

where  C  is  a  normalizing  constant 

is  the  desired  recursive  relation  for  W^.  The  required  density  p  is 
approximately  normally  distributed  and  has  distribution 

P { z ( k+ 1 ) | u ( k)  =  u(i),  Z(k)}  -v  N  {mi(k+l),  C  (k+1)},  (3.2.14) 

where 

m  (k+1)  =  H[*(k)  i(i)(k)  +  r(k)  u(l)(k)]  (3.2.15) 


and 

Ci(k+1)  =  (HM(k+l)  HT  +  R] 

Consider  the  measurement  density  conditioned  on  the  ith  mean 
value  as  given  in  (3.2.14).  This  density  has  covariance  given  by  (3.2.15). 
What  characterizes  the  different  target  "states"  is  the  set  of  N 
Gaussian  curves  used  to  model  the  switching  input.  However,  the  target 
dynamics  remain  the  same  for  all  the  "states".  Consequently,  if  the 


35 


Figure  1.2.2  Typical  Channel  Block  Diagram  of  the  Adaptive  Filter 


process  and  measurement  noise  covariances  Q(k+1)  and  R(k+1) ,  respec¬ 
tively,  are  assumed  to  remain  constant  as  the  target  switches  from 
one  "state"  to  another,  then  none  of  the  quantities  on  the  right  is 
conditioned  on  i.  Under  these  conditions,  for  a  given  value  of 

(k+1),  C^k+l)  has  the  same  value  for  all  values  of  i  =  1,  2 . 

n.  Indeed  it  is  clear  that  for  each  value  of  (k+1) ,  the  entire 
covariance  analysis  is  identical  for  each  filter  in  the  previously 
mentioned  filter  "bank".  Therefore,  the  bank  of  filters  computationaly 
is  not  much  more  than  that  of  a  single  filter. 

3.5  The  Adaptive  Extended  Kalman  Filter 

In  the  one  step  ahead  prediction  by  the  adaptive  filter,  N 
predicted  estimates  X(k+l| k) ^  are  produced,  each  one  conditioned  on 
a  different  vector  (3.2.9).  In  general,  the  weighted  estimate 

N  -  m 

X(k+1 1 k)  =  l  X(k+l|k)U;  W  (k)  (3.3.1) 

i-1  1 

is  closer  to  the  true  state  vector  than  the  individual  conditioned 
estimates  where  the  weights  W^(k)  are  those  which  are  computed  at  the 
end  of  the  kth  iteration.  By  computing  the  Gain  matrix  (3.2.3)  using 
an  H  matrix  linearized  about (3.3.1)  rather  than  N  such  Gain  matrices 
linearized  about  the  individual  estimates  X(k+l|k)  a  more  accurate 

value  is  obtained.  Since  only  one  H  matrix  (and  one  Gain  matrix)  is  now 
being  used,  the  entire  covariance  analysis  becomes  the  same  for  each 
mean  value  and  a  single  filter  suffices  for  this  portion  of  the  algorithm. 
This  is  how  the  conflicting  requirements  of  the  adaptive  and  Extended 
Kalman  filter  algorithms  are  reconciled.  In  summary,  Equations  (3.2.2), 
(3.2.3),  and  (3.2.4)  describe  the  covariance  portion  of  the  Extended 
Kalman  filter  provided  linearization  takes  place  about  X(k+l|k)  as 
defined  by  (3.3.1). 
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B 


The  state  predict  and  state  update  equations  are  now  given 
with  reference  to  the  nonlinear  measurements. 

-  h^X) 

t2  -  h2(X) 

where  h^,  h2  are  defined  in  (3.1)  and  (3.2),  respectively  and, 
X(k+l|k)(i)  -  $X(kjk)(i)  +  Tu(i). 

The  ith  filter  conditional  measurement  residual  is  defined  as 


2(i) (k+1) 


l'1’  (k+1) 

z  (k+1)  -  h  (X(k+l!k)(i)> 
T1 

i*11 (k+1) 

* 

(i) 

z  (k+1)  -  h2(X(k+l|k)  } 

(3.3.2) 


where 

z  (k+1)  and  z  (k+1)  are  the  noisy  measurements  of  r  and  t_, 

T1  T2 

respectively,  at  time  (k+1).  N  such  conditional  measurement  residuals 
are  computed  corresponding  to  the  N  column  vectors  given  by  (3.2.1). 
Then  the  ith  updated  state  estimate  is  computed  as  follows: 

X(i)(k+l|k+l)  -  X(l)(k+l[k)  -  K{X(k+l|k)}[2(i)(k+l)]  (3.3.3) 

for  each  value  of  i  *  1,  2,  ...,  N. 

In  order  to  complete  the  adaptive  Extended  Kalman  filter 
algorithm,  the  linearized  H  matrix  is  now  developed.  Referring  to 


a (k+1) 


H  has 

the 

form 

3Ti 

30 

0 

0 

3T1 

3z 

3t2 

3t2 

3o 

0 

0 

3z 

0  0 


0  0 


(3.3.4) 


X(k+1 ; k) 


Note  that  z  is  the  relative  vertical  separation,  i.e.  zgQ  in  Figure 
3.1.1.  From  (3.1)  and  (3.2)  we  get 
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1 


where 


(z 

[  - 


+  2H02> 


(3.3.5) 


D 

D 

D 


1 

2 

3 


(P2  +  z2  -  4HQ12  -  4HQiz)1/2 
(p2  +  z2)1/2 


(P2  +  Z2  +  4HQ22  +  4HQ2z)1/2 


Note  chat  Equations  (3.3.5)  are  evaluated  using  the  one  step  ahead 
predicted  states;  implicit  in  this,  too,  is  that  the  predicted  values 
for  and  Hq2  are  also  calculated  and  used  in  (3.3.5).  Equations 
(3.2.2),  (3.2.3),  (3.2.4)  and  (3.3.1)  -  (3.3.3)  constitute  the 
adaptive  Extended  Kalman  filter  algorithm. 

A  weighted  estimate  of  the  updated  conditional  estimates  is 
next  formed  using  the  weights  given  by  (3.2.11).  It  should  be  pointed 
out  that  because  the  measurements  now  form  a  vector,  namely 


Z(k+1) 


T^k+1)  +  V1(k*l) 

T,(k+1)  +  v? (k+1) 


when  r.::3  additive  zero  mean  white  Gaussian  measurement  noise, 

the  multivariate  Gaussian  distribution  must  be  used  to  compute  the 
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above  weights.  This  distribution  has  the  mean  value  vector 


“h1(X(k+l|k)(1))“ 

_h2(X(k+l|k)(i)). 

and  covariance  matrix 

C(k+1)  -  [H{X(k+l | k) }P(k+l | k)H{X(k+l j k) }T  +  R(k+1)] 

This  completes  the  discussion  of  the  Adaptive  Extended  Kalman  filtering 
algorithm.  In  the  next  section  specific  design  values  will  be 
discussed  and  the  resulting  filter  performance  evaluated. 


3.6  Specific  Filter  Designs  and  the  Resulting  Filter  Performance 

Using  Synthetic  Data 

The  numerical  values  used  for  submarine  velocities,  depths,  etc 
in  Che  remainder  of  this  report  should  not  be  taken  as  being  representa¬ 
tive  of  actual  values  for  modern  submarines,  but  are  used  exclusively 
for  reference  purposes. 

The  parameter  selection  process  concerns  itself  with  the 
following  parameter  set 


V  »•  Vmax’  *■  V 


where 


a  is  Che  assumed  drag  coefficient 

cjc  is  the  standard  deviation  of  the  correlated  process 

N  is  the  number  of  levels  (mean  values) 

Vmax  iS  the  aS3ume<*  maximum  possible  speed  of  the  target 
set  to  be  tracked 


where  is  the  correlation  time  constant  of  the  correlated  process, 
and  »  £(u-  (♦u-u<V  is  the  mismatch  term  in  the  covariance 

calculation. 

The  set  of  parameter  values  is  summarized  below: 


T  ■  10.0  secs 

a  •  0.05 

a  *  0.18 

c 

0 

°c  -  0.03 

z 

N  -  7 

V  “  ±24  ft/sec 

max 

a  -  0.25 
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-1.2,  -0.8,  -0.4,  0.0,  0.4,  0.8,  1.2 


0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0 


0.06  0 


The  following  initial  conditions  are  used: 
Hq2  -  2000' 


H  »  6000* 

V 


With  a  water  depth  Hw  of  6000',  these  initial  conditions  place 
the  Observer  and  Source  at  depths  (below  the  water  surface)  of  1000'  and 
600',  respectively,  with  the  Source  400'  above  the  Observer.  These 
depths  are  more  realistic  chan  those  used  before.  Because  of  tempera¬ 
ture  gradients  and  a  variety  of  other  reasons,  the  velocity  of  sound 
in  water  is  not  a  constant  but  varies  with  position.  In  an  attempt 
to  take  this  into  account,  the  velocity  is  modeled  as  a  Gaussian  random 
process  with  a  mean  value  C  *  4950  (ft/sec).  This  is  the  modeling 
technique  which  is  used  in  [Si.  The  above  conditions ,  which  provide 
for  a  much  more  realistic  environment  in  which  to  test  the  filter, 
yield  the  following  range  of  time  delays: 


26  a. sec.  -  11  a. sec. 

x 2 .  951  a. sec.  -  471  a. sec. 

With  additive  aeasureaenc  noise  of  5  a. sec.  standard  deviation,  the 
SNR  for  the  aeasureaencs  gets  lower  as  time  progresses  due  to  the 
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decreasing  value  of  t^.  This  decrease  in  (and  to  a  lesser  degree 
t^)  is  caused  by  the  decrease  in  the  difference  between  the  lengths 
of  the  direct  and  surface  reflected  paths  as  the  separation  of  Observer 
and  Source  increases.  In  addition,  since  these  time  delays  are  gen¬ 
erated  using  a  random  velocity  of  sound  in  water,  the  SNR  of  is 
effectively  decreased  even  further. 

Figure  3.6.1  shows  the  relative  positions  of  the  Source  and 
Observer.  Vg  is  the  Source  velocity  with  respec  .  to  the  fixed  X  axis 
and  Vq  is  the  Observer  velocity,  also  with  respect  to  the  fixed  X  axis. 

In  Figures  3.6.2  -  3.6.4,  the  Source  has  a  horizontal  velocity 

V  of  20  (ft/sec)  and  the  Observer  a  horizontal  velocity  V  of  4  (ft/ 
s  o 

sec).  This  yields  a  16  (ft/sec)  relative  horizontal  velocity  of  the 
Source  with  respect  to  the  moving  Observer^  at  a  time  unknown  to  the 
Observer  the  target  makes  a  major  speed  change  as  shown  in  Figure 
3.6.3.  The  percent  error  in  radial  position  is  generally  well  within 
±3%  in  Figure  3.6.2.  The  radial  velocity  estimates  in  Figure  3.6.3 
continue  to  oscillate,  but  with  smaller  excursions  than  with  N  *  3 
levels.  In  Figure  3.6.4,  the  excursions  in  the  z  estimates  are  on 
the  average  in  error  by  50  feet. 

It  appears,  therefore,  that  increasing  the  number  of  mean  values 
from  3  to  7  just  about  maintains  the  same  filter  performance  in  the 
face  of  reduced  SNR  on  and  the  modeling  of  the  velocity  of  sound 
in  water  as  a  random  process  to  generate  the  sonar  time  delay  data. 

The  oscillations  in  radial  rate  p  are  still  present  and  a  slight 
decrease  has  occurred  in  the  quality  of  the  vertical  separation  estimate. 
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Relative  Positions  of  Source  and  Observer 


Estimated  Source  Relative  Radial  Velocity  (3  levels) 


fee 


Figure  3.6.4  Estimate  of  Source  Relative  Vertical  Position 


3.7  Conclusion 


Depending  on  Che  horizontal  distance  separating  the  Observer 
and  Source,  the  ray  paths  in  Figure  3.1.1  may  undergo  multiple  reflec¬ 
tions  before  arriving  at  the  Source.  By  this  is  meant  that  the  surface- 
reflected  wave  may  be  reflected  back  to  the  ocean  floor  and  then  back 
to  the  surface  before  arriving  at  the  Observer.  A  similar  situation 
exists  for  the  bottom-reflected  wave.  Another  problem  which  the  sonar 
time  delay  tracker  may  have  to  contend  with  is  the  problem  of  inter¬ 
mittently  available  time  delays.  For  example,  the  analysis  in  this 
chapter  presupposes  that  the  two  time  delays  are  available  simulta¬ 
neously.  Under  certain  conditions,  this  fortuitous  situation  may  not 
occur  and  the  time  delays  may  be  individually  available  only  at  random 
times,  or  one  or  both  of  them  may  not  be  available  at  all  for  certain 
periods. 

Thus  the  sonar  time  delay  tracker  is  susceptible  to  a  variety 
of  debilitating  influences.  The  decoupling  of  the  bearing  tracker  from 
the  processing  of  the  time  delays,  which  the  Adaptive  Polar  Kalman 
filter  developed  in  this  chapter  permits,  enables  bearing  measurements 
to  be  processed  even  if  the  sonar  time  delays  develop  poor  quality 
or  are  interrupted.  With  the  rectangular  filters,  the  processing  of 
the  bearing  measurements  is  coupled  to  the  processing  of  the  sonar  time 
delays,  and  any  interruption  of  the  latter  brings  the  bearing  estimation 
to  a  halt  as  well.  Thus  filtering  in  polar  coordinates  not  only  reduces 
the  order  of  the  measurement  linearization  matrix,  but  also  localizes 
to  the  p-z  plane  the  deleterious  effects  of  poor  quality  sonar  time 
delay  measurement  data. 
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The  decoupling  of  the  bearing  channel  as  a  result  of  filtering 
in  polar  coordinates  has  two  important  results:  namely,  that  the  order 


of  the  linearized  measurement  matrix  H  is  reduced  by  one  and  the 
effects  of  poor  or  intermittently  available  sonar  time  delays  are 
confined  to  the  p-z  plane.  The  estimates  of  the  Source  radial  velocity 
become  very  poor  under  low  SNR  conditions  on  the  sonar  time  delays. 
Therefore  a  search  for  an  Independent  method  of  estimating  the  Source 
radial  velocity  is  worthwhile.  This  method  should  function  reliably 
and  accurately  regardless  of  the  quality  of  the  sonar  time  delay 
measurements.  The  remaining  part  of  this  report  is  devoted  to 
developing  an  independent  method  which  precisely  satisfies  these 
requirements. 


Chapter  A 


PASSIVE  TRACKING  OF  SOURCE  RELATIVE  RADIAL  VELOCITY  PROFILE  USING 
THE  DOPPLER  EFFECT  AND  KNOWN  CONSTANT  CENTER  FREQUENCY 

4.1  Introduction 

It  is  shown  at  the  end  of  Chapter  3  that  passive  tracking 
using  sonar  time  delays  yields  acceptable  results  only  when  these 
measurements  are  characterized  by  high  signal-to-noise  ratios.  How¬ 
ever,  even  under  such  favorable  conditions,  oscillations  occur  in 
the  estimates  of  the  Source  relative  radial  velocity;  these  oscillations 
increase  as  the  SNR  ratio  decreases.  This,  in  turn,  causes  the  per¬ 
cent  error  in  the  Source  relative  radial  position  to  increase  and  also 
leads  to  much  poorer  estimates  of  the  Source  relative  vertical  position. 
While  it  is  true  that  the  performance  of  any  filter  will  suffer  under 
adverse  SNR  conditions,  the  effect  is  particularly  pronounced  with 
both  the  Extended  and  Adaptive  Extended  Kalman  filters,  because  of 
the  approximate  nature  of  these  filters. 

There  exists  another  and  to  a  considerable  extent  an  indepen¬ 
dent  method  of  determining  estimates  of  the  Source  radial  velocity 
(note  that  the  term  "relative"  is  dropped;  from  this  point  onward, 
and  unless  stated  otherwise,  whenever  Source  position  or  velocity 
is  mentioned,  it  is  implied  that  they  are  relative  to  the  Observer). 

This  method  uses  what  is  commonly  known  as  the  Doppler  effect.  The 
next  section  serves  as  a  brief  introduction  to  the  Doppler  effect. 

The  remainder  of  this  chapter  is  then  devoted  to  using  this 


50 


T 


effect  (assuming  the  "center  frequency"  fQ  is  known)  to  passively  track 
the  radial  velocity  profile  of  a  maneuvering  Source.  The  next  chap¬ 
ter  will  address  this  problem  assuming  the  center  frequency  is  not 
known,  hut  rather  a  random  process. 

4.2  The  Doppler  Effect  for  Sound  Waves 

Consider  Figure  4.2.1(a).  In  this  figure  the  stationary  sound 
source  A  is  transmitting  a  sound  wave  of  frequency  f  (also  called  the 
"center  frequency").  Because  of  his  motion,  the  frequency  of  the 
sound  heard  by  the  listener  B  is  not  the  center  frequency  f  but  some 
other  frequency  £  to  be  determined.  Before  proceeding  further,  it 
should  first  be  pointed  out  that  r  in  Figure  4.2.1  refers  to  the 
separation  of  the  source  and  listener  in  three-dimensional  space — the 
spherical  radius.  For  example,  in  (a)  the  listener  might  be  a  subma¬ 
rine  and  the  source  a  floating  sonar  buoy  providing  navigation  infor¬ 
mation  by  means  of  the  Doppler  effect.  Thus  r  need  not  be  horizontal 
as  drawn.  In  addition,  r  is  the  component  of  the  listener's  velocity 
vector  acting  along  the  radial  direction  r — the  spherical  radius 
race.  With  the  meaning  of  r  and  r  firmly  established,  refer  again 
to  Figure  4.2.1(a),  where  the  listener  is  moving  away  from  the 
stationary  source  with  a  radial  velocity  r.  This  radial  velocity 
imparts  a  shift  [12] 

Af  *  -  f  r/C  (4.2.1) 

o 

to  the  center  frequency,  where  C  is  the  velocity  of  sound  in  the 
surrounding  medium.  This  frequency  shift  is  then  added  to  the  center 
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Figure  4.2.1(b)  Moving  Source  A  and  Stationary  Listener  B 


frequency  and  it  is  the  sum  of  these  two  frequencies  which  the  listener 
hears 

f  r 

ft  *  f0  +  (-  -§“>  -  fQ  [1  -  r/C]  (4.2.2) 

Equation  (4.2.1)  implies  that  for  a  positive  (negative)  r,  the  received 

frequency  f  in  (4.2.2)  is  less  (greater)  than  the  center  frequency  f  . 

t  o 

Now  consider  Figure  4.2.1(b)  where  the  source  is  moving  and 
the  listener  is  stationary.  The  exact  expression  for  the  shift  in  this 
case  is  different  from  that  given  by  (4.2.1).  However,  if  the  speed 
of  the  source  is  small  compared  to  the  speed  of  sound  in  the  surround¬ 
ing  medium,  then  the  shift  imparted  to  f  in  case  (b)  simplifies  to 
that  given  in  (4.2.1)  [12].  This  implies  that  f  in  case  (b)  is  also 
given  by  (  4.  2.  2)  . 

Now  the  typical  underwater  scenario  where  both  the  Source  and 
Observer  (Listener)  are  moving  does  not  fall  into  either  case  (a)  or 
(b)  in  Figure  4.2.1.  To  determine  the  formulas  for  the  Doppler  shift 
in  this  situation,  Equation  (20-11)  in  reference  [12]  states  that  for 
a  moving  Observer  and  Source 


f 


t 


(1  +T) 

where  "the  upper  signs  correspond  to  the  Source  and  Observer  moving 
along  the  line  joining  the  two  in  the  direction  toward  the  other, 
and  the  lower  signs  in  the  direction  cajay  from  the  other."  These 
two  cases  are  now  considered  individually. 
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lull  IUIW 


Case  I  Source  and  Observer  moving  toward  each  other  with  speeds 
vg  and  v  ,  respectively. 


(1  +  -§■)  v  v 

ft  *  fo  - 7--fo<1+f)<1+“#) 

(1  — T-) 


if  -—  <<  1*  Therefore 


ft  =  fo (1  +  t  +  -  f0  (1  +  V-SL-rJL)  (4-2-3) 

Let  r  be  defined  as  the  rate  of  increase  in  the  distance  r  separating 

the  Observer  and  Source.  Then  since  v  and  v  are  positive  numbers 

o  s 

(speeds) 

*  "  "  (vo  +  V 

is  the  relative  radial  velocity  of  the  Source  with  respect  to  the 
moving  Observer  and  is  negative.  Substituting  this  into  (4.2.3) 

ft  *  f0(l  -  r/C)  (4.2.4) 

for  r  <  0. 


Case  IX 


Source  and  Observer  moving  away  from  each  other  with  speeds 
vg  and  v  ,  respectively. 

v 


f 

o 


(1+T> 


(1 


Therefore 


ft  a  fo  (1 


fo  (1 


V  +  V 


(4.2.5) 
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Defining  r  as  Che  race  of  increase  in  r,  chen 

r  -  (v  +  v  ) 

0  s 

and  r  is,  once  again,  Che  relaCive  radial  velociCy  of  Che  Source  with 
respecC  Co  Che  moving  Observer  and  is  now  positive.  SubsciCuCing  chis 
inco  (4.2.5) 

f£  *  fQ  (1  -  r/C)  (4.2.6) 

for  r  >  0. 

Equacions  (4.2.4)  and  (4.2.6)  are  valid  for  r  <  0  and  r  >  0, 
respecCively .  However,  Chese  two  equacions  are  idencical  and  Che 
conclusion  is  chac  eicher  of  chese  is  Che  expression  for  che  received 
frequency  as  measured  ac  Che  moving  Observer  when  r  is  interpreted  as 

the  relative  radial  velocity  of  the  moving  Source  with  respect  to  the 

moving  Observer.  Therefore  Che  sicuacions  in  Figure  4.2.1  which  lead 
Co  Equacion  (4-2.2)  are  boch  special  cases  of  che  more  general  case 

involving  mocion  by  boch  che  Observer  and  Source.  In  summary)  che 

Doppler  shifc  Af  is  given  by  (4.2.1)  where  r  is  incerpreced  as  che 
relacive  radial  velocicy  of  che  Source  wich  respecC  Co  che  moving 
Observer  and  che  received  frequency  measured  ac  che  Observer  is 
given  by  Equacion  (4.2.6).  Having  discussed  che  Doppler  effecC  in 
Che  presence  of  a  moving  Observer  and  Source,  che  use  of  this  effecc 
in  passive  cracking  is  now  discussed. 


4.3  Velocicy  Tracking  Using  the  Doppler  EffecC 

Consider  the  Observer-Source  situation  shown  in  Figure  3.6*1. 
The  distance  r  is  given  by 


1/2 


r  *  <P  +  2so  > 
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and 


*-*<£>  +  %o(-f>  (4-3*1) 

For  the  reasons  cited  in  Chapter  3,  there  are  only  relatively  short 

periods  when  z  /  0.  Also  for  large  values  of  r, 
so 


and 


c-f  >  *  0 


-i 


(4.3.2) 


The  approximations  in  (4.3.2)  have  good  physical  justifications.  The 
Source  and  Observer  are,  because  of  pressure  buildup,  very  restricted 
in  how  far  below  the  ocean  surface  they  can  safely  penetrate.  This 
restriction  on  vertical  movement  in  turn  places  limitations  on  the 
extent  of  their  vertical  separation,  zgQ.  There  is  no  such  restriction, 
however,  on  the  horizontal  distance  separating  them.  Therefore,  while 
the  vertical  separation  z  is  measured  in,  say,  hundreds  of  feet, 
the  horizontal  distance  P  is  measured  in  miles  or  tens  of  miles. 

In  view  of  Equation  (4.3.2),  Equation  (4.3.1)  can  be  approximated  as 
r  -  p  (4.3.3) 

and 

-f 


if  s  - 


fo  P 


, _ o.  • 

(  c  ? 


(4.3.4) 


where  C  is  the  velocity  of  sound  in  water. 

This  simplification  of  Equation  (4.3.1)  to  (4.3.3)  has, 
through  Equation  (4.3.4),  shown  that  the  Doppler  shift  is  a  linear 
function  of  the  radial  rate  p.  Consequently  the  radial  channel  model 
for  the  maneuvering  Observer-Source  dynamics  developed  in  Chapter  2 
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can  now  be  used  Co  passively  Crack  Che  Source  radial  velocicy  using 
(4*3.4).  This  model  is  repeaced  below  for  convenience. 


(A-Tf 

- 
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so 
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(E-l) 
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o  so 
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Equacion  (4.3.4)  shows  chac  Che  radial  velocicy  informacion  is 
concained  in  Che  Doppler  shift.  However,  Che  measured  frequency  f 

“k+l 

is  noc  che  shifc,  buc  racher  Che  sum  of  Che  shifc,  cencer  frequency 
f  and  measuremenc  noise: 


"k+l 


f  +  Af.  ,,  +  v, 
o  k+1  k+1 


(4.3.5) 


where 


Af 


k+1 


-f 

/  O',  ' 

(  C5  Pk+1 


(4.3.6) 


If  che  known  cencer  frequency  fQ  appearing  in  (4.3.5)  is  looked  upon 

as  a  fixed  measuremenc  bias  added  onco  the  shifc  Afj^,  Chen  che  bias 

can  be  direcCly  removed  since  f  is  known  in  Chis  chaDCer.  Therefore. 

o 

by  prefiltering  Che  measured  freauencv  f  ,  the  resulting  prefiltered 

"k+l 

measurement  Af.  . ,  +  v. 

k+1  k+1 
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is  a  linear  function  of  the  state  variable  p  as  given  by  (4.3.6). 


Zk+1  "  “k+l  +  vk+l 


(4.3.7) 


and  the  adaptive  Kalman  filter  algorithm  can  be  applied  to  these  noisy 
measurements  (4.3.7)  using  the  earlier  parameter  values  of  Chapter  3. 

Consider  now  Figure  4.3.1  where  the  Source  and  Observer  are 
moving  horizontally,  the  latter  at  a  constant  velocity  VQ  -  2  ft/sec. 
The  Source  makes  several  abrupt  changes  in  its  velocity  as  shown  by  the 
solid  line  in  Figure  4.3.2.  The  dotted  line  in  this  figure  gives  the 
weighted  estimate  produced  by  the  adaptive  Kalman  filter  using  the 
noisy  frequency  measurements  given  by  (4.3.5).  This  weighted  estimate 
is  seen  to  closely  track  the  abruptly  changing  velocity  profile  of 
the  Source  with  relatively  small  lag  time.  These  good  velocity 
estimates  provide  a  reassuring  test  of  the  accuracy  of  the  state 
equation  for  p  in  compensating  for  Observer  motion.  To  appreciate 
this  point,  first  note  that  p  is  the  relative  radial  velocity  of  the 
Source  with  respect  to  the  moving  Observer.  Consequently,  any  sig¬ 
nificant  error  in  compensating  for  the  component  of  the  Observer  s 
own  velocity  acting  in  the  radial  direction  will  ultimately  result 
in  a  bias  on  the  estimate  of  the  Source  radial  velocity.  This  bias 
i3  clearly  absent  in  Figure  4.3.2. 


Figure  4.3.1  Relative  Position  of  Observer  and  Source  in  Obtaining 
the  Results  in  Figures  4.3.2  and  4.3.3 
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Figure  4.3.2  Relative  Radial  Velocity  of  Maneuvering  Source  with  Respect  to  the  Moving 
Using  Doppler  Frequency  Measurements  and  Known  Center  Frequency 
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The  trajectory  in  Figure  4.3.2  also  provides  a  check  on  the 
"integrating  effect"  of  the  state  model  in  the  following  manner.  A 


very  good  initial  estimate  of  the  Source  radial  position  is  assumed  at 
t  *  0  in  Figure  4.3.3.  Then  the  state  equation  for  "integrates" 

the  combined  effect  of  the  Source  radial  velocity  estimate  and  the  Ob¬ 
server  radial  velocity  to  arrive  at  the  Source  radial  position  estimate 
p  at  each  time  step.  The  percent  error  in  Source  radial  position  at  each 
step  is  then  calculated  and  the  results  are  plotted  in  Figure  4.3.3.  It 
is  seen  that  over  a  period  of  50  minutes  the  percent  error  in  Source 
radial  position  never  exceeds  ±2%.  These  results  attest  to  the  accuracy 
of  the  state  equation  for  p ,  particularly  in  its  compensation  for  Ob¬ 
server  motion. 

4.4  Conclusion 

The  tracking  of  a  maneuvering  Source  velocity  profile  using  the 
Doppler  effect  and  known  center  frequency  fQ  is  a  relatively  straight¬ 
forward  filtering  problem  in  polar  coordinates.  The  adaptive  polar 
filter  developed  in  Chapter  2  performs  well  in  responding  to  abrupt 
changes  in  Source  radial  velocity.  Assuming  a  good  initial  estimate  of 
Source  radial  position,  the  adaptive  filter  "integrates  out"  the  velocity 
estimates  to  yield  good  position  estimates.  These  radial  position  esti¬ 
mates  remain  within  ±2%  of  the  true  radial  position  in  the  presence  of 
several  velocity  changes  by  the  Source.  These  results  indirectly  attest 
to  the  accuracy  of  the  linearized  polar  state  variable  model  for  the 
maneuvering  Observer-Source  scenario,  particularly  in  its  compensation 
for  Observer  velocity. 
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Chapter  5 


PASSIVE  TRACKING  OF  SOURCE  RADIAL  VELOCITY  USING  THE  DOPPLER  EFFECT 
IN  THE  PRESENCE  OF  A  RANDOMLY  VARYING  CENTER  FREQUENCY 

5.1  Introduction 

In  Chapter  4,  the  Doppler  effect  is  used  to  passively  track  the 
Source  radial  velocity  assuming  the  center  frequency  f  is  known.  Typi¬ 
cally,  a  submarine  emits  a  broad  sound  spectrum  having  several  discrete 
line  spectra.  These  line  spectra  are  often  generated  by  machinery  on 
board  the  Source.  As  previously  discussed  any  of  these  frequencies  can 
serve  as  the  center  frequency  fQ  emitted  by  the  Source.  However,  under 
realistic  tactical  situations,  this  center  frequency  fQ  is  either  unknown 
or  a  slowly  varying  random  process.  From  previous  work  we  know  that  the 
Doppler  information  is  contained  in  the  bki&t  Af.  Unfortunately  it  is  not 
the  shift,  but  rather  the  sum  of  the  shift  and  center  frequency  which  is 
measured  by  the  Observer,  as  indicated  in  Chapter  4.  Therefore  the 

measured  frequency  f  is  virtually  useless  for  estimating  the  Source 

m 

radial  velocity  unless  the  center  frequency  is  known  or  estimated.  The 
received  frequency  f^  (ignoring  measurement  noise)  as  measured  by  the 
Observer  is  given  by 

f  p 

f  -  f  +  Af  -  f  -  ~  (5.1.1) 

t  o  o  C 

where  the  parameters  fQ,  p  and  C  have  been  described  in  detail  in  the 
previous  chapter.  Since  fQ  is  unknown,  then  for  any  given  value  for  ft, 
there  exists  an  infinity  of  possible  combinations  of  values  for  f  and  Af 
which  sum  together  to  give  ft>  Any  algorithm,  therefore,  designed  to 


process  f  must  be  cognizant  of  this  fact  and  take  steps  to  select  the 
correct  pair  of  values.  Another  important  point  to  note  is  that  in  order 
to  estimate  fQ,  it  must  be  defined  as  a  state  variable.  The  full  impli¬ 
cation  of  this  can  be  realized  by  referring  to 


Thts  aquation  thorn  that  the.  shl&t  Af  now  involves  the  product  o&  the 
state  variables  p  and  f  and  consequently  a  nonlinear  estimation  problem 
must  be  solved  in  otiden.  to  process  the  VoppleA  measurement,  in  short,  the 
ssumption  that  f  is  unknown  or  random  has  the  following  two  results: 

1.  It  increases  the  number  of  state  variables 
(and  hence  the  order  of  the  system)  by  one  and 

2.  it  increases  the  complexity  of  the  problem  by 
transforming  the  linear  estimation  problem  of 
Chapter  4  into  a  nonlinear  estimation  problem 
involving  the  product  of  state  variables. 

It  was  stated  in  Section  2.2  that  "the  Extended  Kalman  filter  is 
usually  the  first  filter  to  be  used  when  confronted  with  either  system  or 
measurement  nonlinearities."  Rather  than  using  the  Extended  Kalman  filter, 
an  entirely  different  approach  is  used.  This  new  approach  consists  of 
performing  a  "transformation"  on  the  problem  which  shifts  the  nonlinearity 
in  such  a  manner  that  it  can  be  "disposed  of"  in  a  relatively  easy  manner. 
Then  a  "pseudo-linear"  filter  is  developed  to  process  the  nonlinear  measure¬ 
ments  (5.1.1),  using  essentially  basic  linear  filtering  theory.  The  details 
will  become  clear  as  the  development  proceeds. 

5.2  State  Variable  Model  for  Doppler  Measurements  of  a  Random  Center 
Frequency 

The  received  frequency  ft  neglecting  noise  at  time  t  ,  and 

ix'i  X  lv i  X 
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random  center  frequency  fl  becomes 

°K+1 


-f ' 

z  _  f i  ,  /  °k+lN  • 

t  fo  +  C  ^  ^k+1 

k+1  k+1 


(5.2.1) 


Now,  by  using  equation  (2.3.12)  for  Range  rate  (P^+i)  we  can  rewrite  (5.2.1) 
when  the  substitution 


p.  -  f.  -  f' 
k  t,  o. 


-  f' 


ok/C 


is  made. 


After  making  several  cancellations 


f* 


-  f :  +  (-rr^)  e  (f  -  f: )  +  ( 


-f ' 
o 


k+1  k+1 


ck  A 


k+ls  rn.  ' 

c - }  [Fws  + 

pk 


Letting  the  center  frequency  at  time  (k+l)T  be 

f'  »  f'  +  <$f ' 

°k+l  °k  °k 

where  6f°  is  the  purely  random  variation  which  f^  undergoes  to  produce 
k  k 

f'  .  Note  that  ■5f'  is  entirely  unrelated  to  the  Doppler  effect. 
k+1  k  r 

i  °k+l 

Now  for  a  small  variance  on  5fn  ,  the  ratio  — —  on  the  average 


is  near 


k 


uni  ty . 


Thus  by  making  this  approximation  in  (5.2.2),  and  bringing  f'  over  to 

°k+l 

the  left  side  then  defining 


Afk+1  E  \+l  "  f°k+l 


Af,  =  f  -  f' 
k  k  °k 


we  get  for  the  first  Doppler  state  equation 


-f ' 

lfk+i ' E  4fk +  A 


EV’  +  AUs  +  (E-l)  (V0  Co,6do)k-«»s 
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(5.2.3) 


-f '  -f ' 

°k  °k+l 

Note  that  the  approximation  (— — )  was  substituted  for  ( — - )  on  the 

right  side  in  (5.2.3).  In  order  to  make  the  problem  tractable,  the 

assumption  that  f'  is  Gaussian  with  a  mean  value  f  is  made.  Let  f'  be 
r  o  o  o 

described  by  the  following  autocorrelation  function 

2  2  -af  ^  (5.2.4) 

Rf , (t)  «  fQ  +  of  e  o 
o  o 

where 

f  ■  E{ f ' ( t) } 
o  o 

By  converting  to  discrete  time  kT,  (k+l)T,  ...  the  random  process 

W£  ( t,  . . )  becomes 
f  k+1 
o 


”af  ^  1  -af  T 

w'  =  e  ow'  +  — —  (1  -  e  o  )  wc 

°,  o.  f  o. 

k+l  k  o  k 


(5.2.5) 


and  letting 
equation  5.2.5  becomes 


1  ~af  ^ 

J  =  -  (1  -  e  o  )  as  before 

3f 

o 


-a  T 

f'  *  f  +e  ^o  w'  +  J  w, 
o,  ,i  o  f  f 

k+1  °k  °k 


(5.2.6) 
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The  expression  (5.2.6)  is  not  yet  in  proper  state  variable  form 


because  of  the  presence  of  the  unknown  mean  value  fQ.  This  mean  value  can 


be  eliminated  by  subtracting 

> 

k  \ 

and  substituting  this  expression  into  (5.2.6)  for  f  yielding 


f' 

o 


k+1 


-a  T 

f'  -  w’  +e  fo  wl 

o.  f  f 

k  °k  °k 


+  J  w, 


Collecting  like  terms,  the  following  discrete  time  state  equation  for 
the  randomly  varying  center  frequency  f^  becomes 


f ' 
o 


k+1 


-a  T 

f'  +  (e  o  -  1)  w’ 
o,  f 

k  o. 


+  J  w , 


(5.2.7) 


k  "k 

The  state  equations  (5.2.3),  (5.2.5)  and  (5.2.7)  collectively  form  the 
Doppler  frequency  state  variable  model  (5.2.8)  with  measurement  model 
given  by  (5.2.9),  and  are  summarized  on  the  following  page. 

The  following  observations  are  now  in  order.  The  auto¬ 
correlation  function  (5.2.4)  was  chosen  with  the  following  properties  in 

mind:  1.  It  represents  a  random  process  having  a 

mean  value  f  . 

o 

2.  The  exponential  correlation  is  versatile 
in  that  it  provides  a  tractable  yet  real¬ 
istic  model  for  actual  random  phenomena. 

In  addition,  the  parameter  a^  can  be  adjusted 

o 

to  model  processes  having  bandwidths  varying 
from  wide  to  very  narrow. 

The  added  complexity  of  a  random  center  frequency  f'  has  in¬ 
creased  the  order  of  the  state  variable  model  by  one.  Note  that  the  state 
equation  for  f^  (second  row  in  (5.2.8)  has  the  term  J  in  the  ¥  matrix. 

This  term  provides  a  means  of  setting  a  lower  bound  on  the  magnitude  of 

*  (i) 

the  error  covariance  for  the  conditional  estimate  f  .  No  such  term 

o 

exists  in  the  unknown  constant  center  frequency  case.  This  is  an  added 
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"benefit"  of  the  random  center  frequency  assumption. 

Like  the  previous  Doppler  frequency  models,  (5.2.8)  is  not  a 

linear  state  variable  model  because  of  the  presence  of  the  state  variable 

f^  in  the  V  matrix.  However,  since  f^  is  a  random  process,  initialization 

using  the  first  noisy  frequency  measurement  can  be  applied  to  (5.2.8)  to 

make  it  a  pseudo-linear  time-varying  state  variable  model. 

The  presence  of  w'  in  (5.2.8)  indicates  that  this  Doppler 
s 

p 

frequency  state  model  also  requires  interaction  with  the  reduced  order 
radial  channel  model  given  in  Chapter  2  by 
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(i) 
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This  will  allow  w'  (tk)  to  be  estimated  and  then  used  in  the  Doppler  filter 
SP  -af  T 

dynamics  matrix  <p  equation  (5.2.8).  In  addition,  (e  o  -  1)  and 

1  ”af  ^ 

- (1  -  ~  o  )  appearing  in  the  Doppler  frequency  equation  (5.2.7)  are 

af0 

functions  only  of  the  correlation  time  constant 

1 


in  the  autocorrelation  function. 


5.3  Performance  Anlysis  of  the  Doppler  Frequency  State  Variable  Model 
for  the  Randomly  Varying  Center  Frequency  Case 

Figure  (5.3.1)  shows  the  structure  of  the  filter  to  be  analyzed 

in  this  chapter.  Each  block  consists  of  a  Kalman  Filter,  either  Doppler 

frequency  to  the  left,  or  polar  range  estimator  to  the  right.  The  inputs 

are  the  N  discrete  levels  U-  ^  and  each  filter  is  now  separately  computed 

9P 

since  the  state  matrices  are  functions  of  (i) .  Thus  N  sets  of  frequency  and 
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radial  channel  estimates  are  produced,  each  set  being  conditioned  on  a 
different  mean  value  us  ^ ,  i  -  1,  2,  N.  The  reason  why  this 

particular  filter  structure  is  chosen  will  become  clear  later  when  the 
block  labeled  "Adaptive  Estimator"  is  discussed.  For  the  present  the 
operation  of  the  filter  without  the  adaptive  portion  is  examined.  This 
will  serve  to  gain  an  appreciation  for  the  filter's  operating  character¬ 
istics  which  will  then  set  the  framework  for  developing  the  adaptive 
section. 

Figure  (5.3.2)  is  a  block  diagram  of  the  model  used  to  generate 

the  noisy  frequency  measurement  data.  The  quantity  r  appearing  there  is 

the  Course  spherical  radius  rate  and  is  essentially  a  random  process  having 

an  unknown  distribution.  This  rate  then  multiplies  the  random  center 

frequency  f^  and  this  product  is  divided  by  the  velocity  of  sound  in  water. 

In  the  first  series  of  results  to  be  discussed  below,  this  velocity  is 

assumed  to  be  a  constant  C  (switch  open  in  Figure  (5.3.2).  The  switch  will 

then  be  closed  for  the  later  results  when  the  velocity  of  sound  in  water  is 

modeled  as  a  white  Gaussian  random  process  having  a  mean  value  C.  The  zero 

mean  white  Gaussian  curves  on  the  extreme  left  and  right  of  the  figure 

represent  the  distributions  of  wfQk  and  Vj^  in  Equations  (5.2.5)  and  (5.2.9), 

respectively.  Thus  the  noisy  measurement  f  in  Figure  (5.3.2)  is  a 

“k+l 

relatively  complex  combination  of  products,  quotients  and  sums  of  random 
processes,  most  of  whose  distributions  (f^  r,  Af,  f£)  are  unknown. 

The  following  results  are  obtained  using  the  scenario  in  Figure 
(3.6.1)  whereby  the  observer  trails  the  source  with  velocity 

V  -  2  ft/sec,  and  V  *  20  ft/sec. 
o  s 

The  filter  parameters  and  mean  values  are  given  in  Chapter  3.  The  following 
additional  parameters  have  the  indicated  values: 
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5.3.2  Block  Diagram  of  the  Model  Used  to  Generate 
the  Noisy  Doppler  Frequency  Measurement  Data 


■  0.3  Hertz 
o 

■  17  seconds 
o 

f  »  500  Hertz  (mean  value  of  Random  Center  Frequency  X) 

C  *  4950  ft/sec 

Frequency  measurement  noise  standard  deviation  is  the  same  as  before. 
This  latter  set  of  values  is  selected  to  produce  a  random  center 
frequency  having  significant  excursions  about  the  mean  value  of  500  Hertz; 
the  correlation  time  constant  being  less  than  two  sampling  intervals  ensures 
that  the  random  variation  in  frequency  is  quite  rapid. 

The  Source  velocity  of  20  ft/sec  corresponds  to  a  Source  control 

input 

u  =  aV  »  (0.05)(20)  =  1.0 
s  s 

P 

This  Source  input  lies  midway  between  the  sixth  and  seventh  filter  mean 

values  u  ^  *  0.8  and  u  ^  =  1.2,  respectively, 
s  s 

p  p 

Figure  (5.3.3)  shows  the  seventh  (f^^),  sixth  (f^^)  and  fifth 
(f^  ^  ),  virtual  center  frequencies  (VCF)  versus  the  actual  random  center 
frequency  f\  The  remaining  four  VCFs,  which  lie  progressively  further 
below  f  in  the  figure,  are  not  shown  for  reasons  of  clarity.  Note  how 
each  VCF  closely  follows  the  random  variations  of  f^  .  Throughout  the  entire 
time  period  of  18  minutes,  the  actual  center  frequency  is  generally  bracketed 
midway  between  the  seventh  and  sixth  VCFs.  This  is  particularly  satisfying 
since  the  Source  input  is  likewise  bracketed  between  the  sixth  and  seventh 
mean  values . 

When  either  the  Observer  or  Source  changes  its  radial  velocity,  an 
immediate  jump  occurs  in  the  measured  frequency  fm.  Consequently  it  is  im- 
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the  Presence  of  an  Observer  Maneuver 


portant  that  the  filter  be  able  to  discriminate  bet  'een  a  jump  in  the 
measured  frequency  caused  by  an  Observer  maneuver  on  the  one  hand  and  a 
Source  maneuver  on  the  other.  To  investigate  this,  the  Observer  speed  Vq 
is  increased  from  2  ft/sec  to  8  ft/sec  beginning  at  t  *  7.5  minutes.  First 
note  that  since  this  is  an  Observer  maneuver  only,  the  Source  input  has  not 
changed  and  continues  to  be  bracketed  between  the  two  mean  values  as  before. 

From  Figure  (5.3.3)  this  Observer  maneuver  4A  not  reflected  in  the 
estimates  and  the  actual  center  frequency  continues  to  be  bracketed  between 
the  sixth  and  seventh  VCFs  as  before.  This  is  to  be  expected  as  this  shows 
that  the  filter  recognized  the  jump  in  measured  frequency  as  having  been 
caused  by  an  Observer  maneuver. 

Consider  now  Figure  (5.3.4)  which  shows  the  three  virtual  Doppler 
shifts  Af^,  Af^  and  Af^  versus  the  actual  Doppler  shift.  Here,  the 

*  /  £  \  A  /  “7\ 

actual  shift  is  bracketed  between  Af'  ;  and  Afv  as  required  by 

u  ^  <  u  <  u  Note  that  the  Observer  maneuver  ti  reflected  in 

3P  SP  SP 

the  virtual  Doppler  shifts.  To  understand  why  this  is  necessary,  first  note 
that  a  sudden  change  in  the  Observer  velocity  Vq  causes  a  sudden  change  in 
the  Source  relative  radial  velocity  p.  This  sudden  change  in  p  in  turn 
causes  a  sudden  change  in  the  actual  Doppler  shift  Af.  However,  since  the 
Source  input  has  not  changed,  the  virtual  shifts  (VS)  must  change  in  such  a 
way  that  Af  remains  bracketed  midway  between  Afv  and  Af'- 

Figures  (5.3.3)  and  (5.3.4)  together  illustrate  how  the  filter 
exploits  its  degree  of  freedom  by  satisfying  the  required  inequalities.  For 

example,  note  how  the  VCFs  in  Figure  (5.3.3)  satisfy 

:,(5)  <  ^ i (6)  <  £,(7) 

lo  o  o 

while  the  VS's  in  Figure  (5.3.4)  satisfy 

Af<5)  aV6)  aV7> 
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Actual  versus  Virtual  Doppler  Shifts  In 
esence  of  an  Observer  Maneuver 


The  actual  Source  radial  velocity  p  in  Figure  (5.3.5)  remains 

bracketed  between  the  virtual  radial  velocities  (VRV)  p^^  and  p^  and  p^^ 

in  the  presence  of  the  Observer  maneuver.  Here  again,  since  the  Source 

input  has  not  changed,  the  filter  responds  by  changing  the  VRVs  to  maintain 

the  actual  value  midway  between  the  sixth  and  seventh  VRVs. 

The  behavior  of  this  filter  in  the  presence  of  a  Source  maneuver 

is  examined  in  the  next  three  figures.  With  the  same  set  of  parameters  and 

initial  conditions  as  in  the  preceding  figures,  the  Source  velocity  V 

changes  from  20  to  12  ft/sec  beginning  at  t  =  7.5  minutes.  The  new  Source 

input  u  =  0.6  is  midway  between  u  ^  and  u  ^ . 

P  p  p 

The  response  to  this  Source  maneuver  is  radically  different  from 

that  for  the  Observer  maneuver.  Consider  for  example  Figure  (5.3.6)  where 
we  see  that  after  the  Source  maneuver  at  t  a  7.5  minutes,  the  VCFs  all 
shift  upward  by  an  amount  sufficient  to  have  the  actual  center  frequency 
f^  bracketed  midway  between  the  fifth  and  sixth  VCFs.  The  response  time  of 
th'-  rilter  is  quite  good  as  this  upward  shift  by  the  estimates  is  immedi¬ 
ately  perceptible  at  the  initiation  of  the  maneuver.  Note  how  the  VCFs 

continue  to  follow  the  random  variations  of  f’  both  before  and  after  the 

o 

maneuver . 


In  Figure  (5.3.7)  the  actual  shift  which  is  initially  bracketed 
between  the  sixth  and  seventh  VS’s  before  the  maneuver  ends  up  midway 
between  the  fifth  and  sixth  VS's  after  the  maneuver.  This  same  situation 
applies  to  the  radial  velocity  p  in  Figure  (5.3.8).  Note  how  the  degree 
of  freedom  continues  to  be  exploited  by  the  filter  in  Figures  (5.3.6)  and 
(5.3.7) . 

To  summarize,  the  filter  bank  in  Figure  (5.3.1)  is  shown  to  be 
capable  of  tracking  the  Doppler  shift,  random  center  frequency  and  radial 
velocity  fluctuations  in  the  presence  of  either  an  Observer  or  a  Source 
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maneuver.  It  has  been  determined  that  the  filter  bank  can  discriminate 
between  a  sudden  change  in  received  frequency  caused  by  a  sudden  maneuver 
on  the  part  of  either  the  Source  or  Observer.  For  an  Observer  maneuver, 
the  estimates  adjust  themselves  to  maintain  the  previously  existing  relative 
positions  of  the  actual  values  versus  the  estimated  values.  For  a  Source 
maneuver,  the  estimates  are  adjusted  to  reflect  the  new  position  of  the 
Source  input  with  respect  to  the  filter  mean  values. 

The  preceding  results  look  at  the  filter's  behavior  from  the  point 
of  view  of  the  individual  conditional  estimates.  Of  course,  in  reality  one 
does  not  know  between  which  of  the  levels  the  actual  Source  input  lies  and 
consequently  an  adaptive  technique  must  be  developed  before  the  filter  can 
hope  to  be  of  any  practical  utility.  However,  in  order  to  develop  this 
adaptive  technique,  the  mean  value  of  the  Source  center  frequently  must 
first  be  known.  The  following  section  discusses  how  to  obtain  this  in¬ 
formation. 


5,4  Determination  of  the  Mean  Value  of  the  Source  Randomly  Varying  Center 
Center  Frequency 

The  determination  of  the  Source  center  frequency  mean  value  is  of 

crucial  importance  to  the  problem  at  hand.  Since  the  VCFs  are  separated  by 

less  than  two  Hertz,  it  is  imperative  that  the  method  chosen  yield  a  very 

accurate  value  for  f  . 

o 

The  technique  consists  of  using  the  Adaptive  Extended  Polar  filter 
developed  in  Chapter  3.  If  the  weighted  estimate  p  produced  by  this  adapt¬ 
ive  state  estimator  is  averaged  to  yield  an  average  value  for  the  steady 
state  radial  velocity,  then  this  information  can  be  used  to  achieve  a  very 
accurate  numerical  value  for  E{f'}.  To  see  how  this  works,  consider 


f 


m 


+  v 
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where  p  is  the  averaged  steady  state  radial  velocity.  Substituting 

L  . 


for  p  and  solving  for  f  we  get 

f 

f _ 5 - 

o  r  (5.4.1) 


where  it  has  been  assumed  in  (5.4.1)  that  v  ■  0  where  v  is  the  average  of 

the  zero  mean  frequency  measurement  noise.  This  method  of  obtaining  f  is 

very  robust,  and  is  illustrated  in  the  next  section. 

5.5  The  Adaptive  Doppler  Frequency-Radial  Channel  Tracking  Filter  for 
the  Randomly  Varying  Center  Frequency  Case 

With  a  method  to  determine  a  numerical  value  f  for  the  Source 

o 

center  frequency  mean  value  now  available,  the  last  remaining  problem  is 

the  development  of  the  block  labeled  Adaptive  Estimator  in  Figure  (5.3.1) 

The  key  to  this  new  adaptive  technique  is  the  implicit  assumption  that 

the  mean  value  otf  the  SouA.ce  xandom  centex  xequency  doeA  not  change  mAh 

time,  ok  vexy  Alowly  ovex  long  pexAodi  o&  time. 

The  following  analysis  presupposes  that  this  mean  value  has  already 

been  determined  by  the  method  presented  in  Section  (5.4.4).  Consider  the 

set  of  conditional  estimates  appearing  at  the  output  of  each  Doppler  frequency 

Kalman  filter  in  Figure  (5.3.1).  If  each  of  the  virtual  center  frequencies 

f'  (1),  i  -  1,  2 . N 

°k 

is  averaged  using  a  sliding  window  of  length  L,  then  the  following  set  of 
average  VCFs  is  obtained: 

f  (i)  s  i  I  f'(1)  ,  i  *  1,  2,  ...,  N  (5.5.1) 

°k  L  y«1  °k-Y+l 


S3 


where  k  is  che  current  iteration  of  the  filter  bank  and  L  is  the  sliding 

window  length.  Now  because  each  of  these  VCFs  faithfully  reproduces  the 

fluctuations  of  f^  during  any  averaging  period,  the  average  value  of  f^ 

is  bracketed  between  some  pair  of  (5.5.1).  This  conclusion  is  simply  a 

necessary  result  of  the  inequalities  if 

u  (i)  <  u  <  u 
p  P  P  then 

J  (j)  <  f  <  }  CD  (5.5.2) 

°k  0  °k 

for  some  set  of  consecutive  integers  (i,  j)  drawn  from  the  set  (1,  2,  . . . ,  N) . 


Mich  the  numerical  value  f  for  the  Source  mean  value  f  in 

o  o 

hand,  it  is  a  relatively  easy  matter  to  determine  the  pair  (i,  j)  after 
each  filter  iteration.  This  information  is  then  used  to  form  the  fol¬ 
lowing  weighted  estimates  of  the  Source  center  frequency,  Doppler  shift 
and  Source  radial  velocity: 


*fk 


f r  <j)  (f  -  7  )  +  f ' 
o,  O,  0  o, 

k  _ k _ k 


) 


Af 


(j) 


(f. 


<i)  _  7 


-  V 


<if£i}  (7 

k  o 


f  «>, 

°k 


V 


+  p. 


(i) 


(f  (i>  -  f  «>) 
°k  °k 


(5.3.3) 


(f.  5.4) 


(5T5.5) 
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Note  that  f'  f'  ...»  p.  ^  are  the  iM>ta.ntaneouA  estimates  produced 

°k  °k  k  _ 

*  *■  '  (£) 
by  the  filter  in  Figure  (5.3.1)  as  distinct  from  the  averaged  estimates  f  . 

k 

etc.  Figure  (5.5.1)  gives  a  block  diagram  of  this  Adaptive  Estimator. 

To  gain  a  further  appreciation  of  how  the  Estimator  works,  let  the  Source 

input  initially  be  bracketed  midway  between  u  ^  and  u  Next  let  the 

s  s 

P  (1)  p 

Source  input  switch  to  a  value  very  close  to  ug  J  .  After  an  initial  trans¬ 
ient,  the  VCF  f'  After  an  initial  transient,  the  VCF  f’  ^  approaches 

°k  °k 

in  value  the  actual  Source  center  frequency  f'  .  Tkii  in  tu/in  meani  that  the 


averaged  l/C F  f  ^  appAoache a  the.  value  f 


which  implies  that 


and 


f  -  f  0 

°  °k 


(i)  J  (J)+(J  (i) 

c  “  °k  °k  *  fo) 


*  (i) 

The  result  is  that  the  coefficient  of  f'  in  the  numerator  of 


(5.5.3)  approaches  zero,  and  the  denominator  approaches  in  value  the  co- 
A  (1) 

efficient  of  f '  in  the  numerator. 

°k 


Therefore  as 


U) 


P  P 

then  the  weighted  estimate 

f »  -►  f  * 

°k  °k  °k  . 

This  same  conclusion  applies  to  the  weighted  shift  and  radial  velocity 
estimates  given  in  (5.5.4)  and  (5.5.5),  respectively.  Note  that  the  Source 


input  u  was  switched  close  to  a  mean  value  merely  for  illustrative  purposes. 
8P 

The  Source  input  can  in  reality  switch  anywhere  within  its  continuum  of 
possible  values  and  the  Adaptive  Estimator  will  respond  appropriately. 
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The  results  in  the  next  series  of  figures  are  obtained  using  the 

scenario  in  Chapter  3.  The  quantities  appearing  have  the  values 

zso  =  feet  vertical  separation 

V  =  4  ft/sec 

o 

Vs  *  20  ft/sec 

u  *  1.0 

s 

P 

The  filter  parameters  are  the  same  as  section  3.6  and  the  following  additional 

parameters  have  the  indicated  values: 

o^  *  0.1  Hertz 
o 

*  17  secs 

o 

f  =  500  Hertz 
o 

C  =  white  Gaussian  random  process  with  a  mean 
value  4950  ft/sec 

Frequency  measurement  noise  standard  deviation  is  the  same  as  before. 

The  switch  in  Figure  (5.3.2)  is  closed  in  order  to  generate  the 
noisy  frequency  measurements  using  a  random  velocity  of  sound  in  water. 
Therefore  this  situation  involving  a  non-zero  vertical  separation  between  the 
Observer  and  Source  is  a  much  more  realistic  scenario  which  when  coupled  with 
the  random  center  frequency  and  random  velocity  of  sound  in  water,  collective¬ 
ly  provide  a  realistic  environment  in  which  to  test  the  Adaptive  Estimator. 

In  Figure  (5.5.3)  the  seventh,  sixth  and  fifth  VCFs  are  shown  together 
with  the  weighted  estimate  produced  by  the  Adaptive  Estimator.  Beginning  at 
t  *  13.5  minutes,  the  Source  velocity  changes  to  Vg  =  14  ft/sec  corresponding 

to  u  *  (0.05)(14)  *  0.7 

s 

P 

which  is  slightly  less  than  ug  ^ .  With  the  Source  input  initially  midway 

between  u  ^  and  u  the  actual  center  frequency  f'  is  seen  to  be 

s  s  o 

P  P 

bracketed  midway  between  the  sixth  and  seventh  VCFs.  At  t  *  13.5  minutes,  the 
VCFs  respond  and  by  t  ■  15  minutes  the  sixth  VCF  has  shifted  up  close  to  f^  . 


To  understand  the  weighted  estimate  plot,  first  remember  that  the 
mean  value  of  the  Source  center  frequency  must  first  be  determined  as  in 
Section  5.4.  This  process  is  taking  place  in  a  brief  time  interval  pre¬ 
ceding  t  -  9  minutes  (but  not  lasting  the  full  9  minutes) .  Consequently 
during  this  period  the  weighted  estimate  is  initialized  to  the  first  measure¬ 
ment,  like  all  the  initial  f^  estimates.  The  quentities  in  Equation  (5.4.1) 
have  the  following  computed  average  values: 

f  =  498.396  Hertz 

m 

p  =  16.134  ft/sec  (obtained  from  averaging 

the  Augmented  Extended  Polar  filter 

estimate) 

Substituting  these  values  into  (5.4.1),  the  numerical  value 
f  »  500.0254  Hertz 

O  T 

is  computed  for  fQ.  Note  that  while  there  is  a  slight  error  in  p  =  16.134 — 

the  actual  Source  relative  radial  velocity  is  16  ft/sec —  the  computed  mean 

value  of  fQ  *  500.0254  Hertz  is  nevertheless  extremely  accurate.  To  see  how 

robust  this  method  is,  consider  for  example  if  the  average  radial  velocity 

estimate  produced  by  the  filter  of  Chapter  3  were  15.0  instead  of  16.134; 

the  computed  mean  value  in  this  case  is 

f  *  499.911 

o 

which  is  still  quite  good.  It  should  be  borne  in  mind  that  the  sonar  time 
delay  measurements  used  by  the  Adaptive  Extended  Polar  Filter  to  produce  p 
are  also  generated  using  a  random  velocity  of  sound  in  water. 

In  Figure  (5.5.3),  the  Adaptive  Estimator  (7.5.3)  is  activated  at 
t  *  9  minutes  with  the  learning  period  completed.  The  weighted  estimate  f’ 

k 

rises  immediately  to  the  actual  center  frequency  and  thereafter  tracks  f^ 
very  closely,  even  in  the  presence  of  the  Source  maneuver  at  t  *  13.5  minutes. 
Note  how  the  weighted  estimate  is  almost  identical  in  value  to  the  actual  f^ 
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after  the  Source  maneuver  is  completed. 

In  Figure  (5.5.4),  the  Adaptive  Estimator  likewise  prodices  a 
weighted  estimate  of  the  Doppler  shift  which  is  practically  identical  in 
value  to  the  actual  Doppler  shift.  The  response  time  of  this  weighted 
estimate  to  the  Source  maneuver  at  t  *  13.5  is  relatively  brief  with  only 
a  minimum  transient  in  the  weighted  estimate  before  learning  the  new  Doppler 
shif  t . 

The  radial  velocity  weighted  estimate  p  in  Figure  (5.5.5)  displays 
all  the  desirable  qualities  of  the  two  preceding  figures.  The  response  time 
to  the  Source  maneuver  is  excellent.  Note  how  p  both  before  and  after  the 
maneuver  is  extremely  close  to  the  true  radial  velocity.  The  vertical  scale 
in  each  figure  provides  a  convenient  gauge  of  the  accuracy  of  the  weighted 
estimates . 

The  Adaptive  Estimator  in  Figure  (5.5.1)  can  also  provide  weighted 

estimates  both  of  the  Source  control  input  u  and  Source  correlated  acceler- 

SP 

ation  term  w'  using  the  same  weighting  technique  applied  to  the  other  state 
P 

variables.  For  example,  in  the  above  figures  the  actual  Source  control  inputs 
before  and  after  the  maneuver  are  1.0  and  0.7  respectively.  The  Adaptive 
Estimator  produces  corresponding  weighted  estimates  having  the  following  range 
of  values: 

Before:  0.988  -  1.07 
After:  0.64  -  0.81 

Note  that  these  radial  velocity,  control  input  and  correlated  acceleration 
weighted  estimates  would  be  of  considerable  utility  in  a  combined  sonar  time 
delay-Doppler  frequency  tracking  filter.  Such  a  combination  would  exploit 
the  best  advantages  of  both  filter  types.  For  example,  with  the  Doppler-radial 
channel  filter  providing  the  above  weighted  estimates,  the  order  of  the  time 
delay  tracker  could  be  reduced  by  two  from  that  in  Chapter  3.  This  reduction 
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Figure  5.5.5  Actual  versus  Weighted  Radial  Velocity  in  the 
Presence  of  a  Source  Maneuver 
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in  order  is  possible  because  the  two  state  variables  p,  w'  appearing  in 

s 

p 

(3.2.1)  are  now  estimated  by  the  Doppler  tracker.  Thus  the  state  vector  for 
the  time  delay  tracker  would  now  consist  of  the  following  four  state 
variables: 


[p 


z 


so 


z 

so 


(5.5.6) 


z 

The  high  quality  estimates  of  p,  w'  and  u  produced  by  the  Doppler  filter 

s  s 

p  p 

should  in  turn  improve  the  performance  of  the  time  delay  tracker.  In 
addition,  by  embedding  the  adaptive  feature  of  the  combined  filter  in  the 
Doppler  section,  the  need  for  the  semi-Markov-based  weights  would  be  elimin¬ 
ated.  This  in  turn  means  that  the  oscillations  in  the  weighted  estimates 
experienced  in  Chapter  3  and  which  are  caused  by  the  weights  switching  back 
and  forth  would  also  be  eliminated.  Note  that  the  simplified  time  delay 
tracker  involves  much  more  than  a  mere  reduction  in  the  size  of  the  state 
vector.  It  also  implies  an  elimination  of  all  but  one  of  the  mean  values 
in  each  channel.  For  example,  the  Source  input  would  not  be  modeled  as  a 
4-tng^e  Gaussian  curve  in  each  of  the  p  and  Z  channels.  The  mean  values  of 
these  curves  would  be  the  Doppler  provided  weighted  estimate  of  the  Source 
input  and  zero,  respectively.  Thus  the  augmented  approach  using  several 
levels  could  also  be  dispensed  with. 

Before  concluding  this  section,  the  following  additional  comments 
need  to  be  made.  If  it  is  known  that  the  Source  center  frequency  is  a  single 
constant  tone,  the  modeling  of  this  center  frequency  as  a  narrowband  process 
is  nevertheless  advised  in  order  to  prevent  the  error  covariance  and  Gains 
from  becoming  too  small  and  causing  divergence.  The  constant  tone  could  be 
accurately  modeled  as  a  highly  correlated  random  process  by  selecting  an 
appropriately  large  (small)  value  for  (a^  ). 


o  o 
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In  obtaining  the  results  discussed  in  this  seccion,  the  pre¬ 
filtered  measurement  for  the  ith  Radial  Channel  Kalman  Filter  in  Figure 
(5.3.1)  is  obtained  by  subtracting  the  ith  VCF  from  the  received  noisy  fre¬ 
quency  measurement. 

Prefiltered  Measurements:  (f  -  f'  ^),  i  ■  1,  2 . N 

°k+l  °k+l 

An  interesting  variation  on  this  would  be  to  use  the  ith  virtual  Doppler 
shift  as  the  prefiltered  measurement 

*  (i) 

Prefiltered  Measurements:  &f  ,  i  -  1,  2,  ...,  N 
This  conceivably  might  improve  the  filter  performance  under  high  measure¬ 
ment  noise  conditions. 

When  the  learning  period  is  over  and  f  has  been  calculated,  a 
considerable  savings  in  computation  can  be  achieved  by  cycling  only  the  two 
(i,  j)  filters  where  the  set  (i,  j)  (was  previously  defined).  When  a 
maneuver  occurs,  the  entire  bank  could  be  activated  until  the  new  set  (i,  j) 
has  been  determined,  at  which  point  (N-2)  filters  could  again  be  deactivated. 
Note  that  detecting  Source  maneuvers  for  activating  the  filter  bank  is  a 
relatively  simple  matter  because  a  Source  maneuver  is  immediately  heralded 
by  a  sudden  shift  in  the  measured  frequency  fn  as  well  as  in  the  estimates 
produced  by  the  two  activated  filters. 

5.6  Conclusion 

This  chapter  presents  possible  the  first  attempt  to  solve  the 
problem  of  passively  tracking  the  radial  velocity  of  a  maneuvering  Source 
using  the  Doppler  effect  in  the  presence  of  a  randomly  varying  center  fre¬ 
quency  . 

This  problem  in  general  is  highly  nonlinear,  particularly  if  one 
is  not  careful  in  selecting  the  polar  coordinate  system  in  which  to  filter. 
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An  Adaptive  Estimator  based  on  the  mean  values  of  the  virtual  center  fre¬ 
quencies  was  developed.  This  Adaptive  Estimator  was  rigorously  tested  under 
realistic  environmental  conditions  which  include  a  rapidly  varying  center 
frequency  and  velocity  of  sound  in  water  modeled  as  a  white  Gaussian  random 
process.  The  structure  of  this  Adaptive  Estimator  is  made  possible  by  con¬ 
ditioning  both  the  Doppler  filter  bank  and  radial  channel  filter  bank 
estimates  on  the  individual  mean  values  u  i  •  1,  2,  ...»  N. 

S 

P 

It  is  believed  that  the  adaptive  filter  developed  in  this  chapter 
represents  the  first  significant  attempt  to  deal  with  the  random  center 
frequency  case  in  all  of  its  complexities.  The  results  obtained  for  a 
scenario  involving  a  maneuvering  Source  are  very  encouraging  and  it  is  felt 
that  the  filter  merits  further  study.  One  feature  worth  looking  at  is  vary¬ 
ing  the  window  length  in  the  Adaptive  Estimator.  The  window  length  used  here 
is  L  -  8,  but  a  shorter  window  length  might  yield  a  still  faster  response. 

Of  course,  a  shorter  window  length  would  also  compound  the  oscillation  pro¬ 
blem  encountered  for  large  o^  ,  so  some  tradeoff  is  called  for.  These  small 

o 

oscillations  could  conceivably  be  eliminated  by  averaging  the  output  of  the 

Adaptive  Estimator  using  a  window  having  a  length  equal  to  a  few  correlation 

time  constants  x^  . 

o 

Another  avenue  of  exploration  is  to  see  if  some  connection  exists 
between  the  mean  value  of  the  center  frequency  and  the  Source  dynamics.  If 
such  a  relationship  actually  exists  in  practice,  then  proper  modeling  dictates 
that  at  least  an  approximate  relationship  be  incorporated  into  the  illter. 

Note  that  in  this  case  the  degree  of  freedom  is  no  longer  at  the  disposal  of 
each  filter,  and  consequently  the  semi-Markov  weights  can  be  used  to  compute 


the  weighted  estimates. 


The  characterisClcs  of  a  combined  Doppler  frequency-sonar  time 


delay  tracking  filter  are  investigated  in  the  next  chapter. 


THE  INTEGRATED  ADAPTIVE  DOPPLER  FREQUENCY- REDUCED 
ORDER  SONAR  TIME  DELAY  TRACKING  FILTER 

6.1  Introduction 

In  Chapter  3  the  Adaptive  Extended  Polar  Kalman  Filter  Is  de¬ 
veloped  to  process  sonar  time  delays.  This  filter  suffers  from  large 
oscillations  in  the  radial  velocity  estimates  even  under  high  signal- 
to-noise  ratios.  These  oscillations  prompted  a  search  for  an  indepen¬ 
dent  method  to  estimate  the  Source  radial  velocity,  culminating  in  the 
Adaptive  Doppler  frequency  tracking  filter  just  presented.  Since  the 
Doppler  measurements  contain  no  information  on  the  Source  radial 
position,  the  Adaptive  filter  in  Chapter  5  can  provide  estimates  only 
of  the  Source  radial  velocity.  Therefore,  a  hybrid  filter  formed  by 
integrating  the  Adaptive  Doppler  and  sonar  time  delay  tracking  filters 
offers  the  chance  to  exploit  the  best  features  of  both  filter  types. 

Thus,  for  example,  the  good  radial  velocity  estimates  produced  by  the 
Adaptive  Doppler  filter  can  be  used  by  the  sonar  time  delay  filter  to 
produce  improved  estimates  of  both  the  Source  radial  position  and 
Source  vertical  separation.  To  combine  the  two  estimators  the  standard 
N  level  Adaptive  time  delay  filter  is  initially  run.  Then  p  is  smoothed, 
and  used  to  compute  f  in  the  Doppler  filter.  Once  this  period  of 
initialization  is  over  the  standard  N  level  filter  is  replaced  by  the 

reduced  order  state  variable  model  where  now  N  >  1  and  u  is  available. 

SP 

The  time  delay  filter  is  derived  in  the  next  section. 


6.2  State  Variable  Model  for  the  Reduced  Order  Sonar  Time  Delay 
Tracking  Filter 

The  Adaptive  Doppler  tracking  filter  developed  in  Chapter  5  provides 
estimates  of  the  following  quantities  from  the  noisy  frequency  measure¬ 
ments  f  : 


p,  w^  and  us 
P  P 


(6.2.1) 


Therefore  the  reduced  order  state  variable  vector  for  the  sonar  time 
delay  tracker  previously  given  is  repeated  here  for  convenience. 


[  P  2  2  W1  ] 

SO  so  s 

z 

Returning  to  the  state  model  (3.2.1)  for  the  time  delay  tracking  fil¬ 
ter,  the  following  state  equations  are  obtained  (the  subscript  sq  on 
z  is  dropped  for  convenience) : 


Vi  '  »k  +  +  K  +CuS  +  <A-T)(Vo  Co.8)  +  D»,  (6.2.2) 

pk  pk  k  pk 


2,  +  Az.  +  Bw'  +  Cu  +  )A-T)  2  +  Dw 


*k+l  *  "k  +  AZk 


(6.2.3) 


zk+l  *  Ezk  +  ^'s  +  Aus  +  (E_1)  Zo.  +  ^s 

zk  zk  k  Zk 


(6.2.4) 


i  -aT  .  ,  T 

w'  *  e  w'  +  J  w 
s  s  s 

\+l  zk  zk 


(6.2.5) 


Since  the  quantities  appearing  in  (6.2.1)  are  estimated  by  the  Doppler 


tracking  filter,  they  can  be  considered  as  "deterministic  inputs"  for 
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Che  reduced  order  sonar  time  delay  filter.  Combining  these  quantities 
with  the  deterministic  inputs  of  (3.2.1)  the  following  expanded  deter¬ 
ministic  input  vector  is  obtained: 


V  Cos  6 
o  so 


(6.2.6) 


Using  this  expanded  input  vector,  the  state  variable  model  (6.2.7)  (next  page) 
emerges  from  the  state  equations  (6.2.2)  -  (6.2.5).  The  corresponding 
linearized  measurement  matrix  H  is  given  in  Equation  6.2.8). 


3T1 

3t 

3p 

Tz 

3t2 

3t 

3p 

Tz 

where  Che  partial  derivatives 

3T1  3t2 

3p  *  3z 


(6.2.8) 


are  defined  in  Equations  (3.3.5  )• 


(•.2.7) 


Note  that  in  the  vector  (6.2.6),  u  is  no  longer  one  of  the 

s 

(i)  P 

preselected  mean  values  ug  ,  i  *  1,  2,  ....  N  but  rather  is  the 

P 

weighted  estimate  of  the  Source  input  produced  by  the  Adaptive  Doppler 

tracking  filter.  In  addition,  and  for  the  reasons  given  in  Section  3.2, 

the  input  ug  has  the  constant  value 
z 

u  *  0. 
s 

z 

Thus  neither  the  vector  (6.2.6)  not  the  state  model  (6.2.7)  is  condi¬ 
tioned  on  i,  and  consequently  only  a  i-ingtn  Extended  Kalman  filter  rather 
than  a  bank  of  N  such  filters  need  be  executed  to  process  the  noisy  time 

delay  measurements  z  and  z  ,  after  the  initialization  or  learning 

T1  T2 

period  is  over. 

6.3  Performance  Analysis  of  the  Integrated  Adaptive  Doppler  Frequency- 
Reduced  Order  Sonar  Time  Delay  Tracking  Filter 

Figure  6.3.1  gives  a  block  diagram  of  the  integrated  Adaptive 

Doppler  frequency-reduced  order  sonar  time  delay  tracking  filter.  In 

this  figure,  the  weighted  estimates  (6.2.1)  produced  by  the  Adaptive 

Doppler  filter  are  fed  to  the  reduced  order  sonar  time  delay  filter 

where  they  are  used  in  processing  the  noisy  time  delay  measurements 

z  and  z 
1  t2 

The  results  presented  in  the  remainder  of  this  chapter  consist 
of  side  by  side  comparisons  between  the  integrated  Adaptive  Doppler  fre¬ 
quency-reduced  order  sonar  time  delay  filter  (hereafter  referred  to  as 
the  hybrid  filter)  and  the  pure  time  delay  Adaptive  Extended  Polar  filter 
developed  in  Chapter  3.  These  results  take  the  form  of  a  series  of  graphs, 
each  graph  containing  superimposed  plots  of  the  correcponding  estimates 
produced  by  the  two  different  filters.  Thus,  a  very  effective  comparison 
can  be  made. 
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Figures  6.3.2  -  6.3.4  are  obtained  using  the  filter  parameters 

of  Chapter  3,  zSq  =  400  feet,  Vo  =  4  ft/sec,  Vs  *  20  ft/sec,  be  *  3  ft/sec 

u  *  1.0,  b,.  =  .lhz,  t,  ■  17  sec  and  fQ  =  500  hz.  In  addition,  the 
ro  ro 

random  velocity  of  sound  in  water  is  used  to  generate  both  the  noisy 
frequency  and  time  delay  measurement  data.  It  should  also  be  pointed  out 
that  the  time  t  *  0  in  these  figures  corresponds  to  the  instant  when  the 
learning  period  of  the  Adaptive  Doppler  frequency  tracking  filter  has 
just  been  completed. 

In  Figure  6.3.2  the  percent  error  in  Source  radial  position  p 
produced  by  the  pure  time  delay  filter  is  larger  by  several  orders  of 
magnitude  than  that  produced  by  the  hybrid  filter.  Indeed,  whereas  the 
hybrid  filter's  errors  are  well  within  ±1%  of  the  true  value,  the  errors 
produced  by  the  other  filter  are  larger  than  4%  for  values  of  t  15 
minutes. 

In  Figure  6.3.3  the  estimates  of  the  Source  radial  velocity  pro¬ 
duced  by  the  hybrid  filter  are  far  superior  to  those  produced  by  the 
pure  time  delay  filter.  While  the  hybrid  estimates  are  generally  within 
±1.5  ft/sec  of  the  actual  value  (16  ft/sec),  the  other  filter's  estimates 
are  in  error  by  as  much  as  ±5  ft/sec.  Note  that  with  the  adaptive  fea¬ 
ture  embedded  in  the  Doppler  tracker,  the  hybrid  filter  estimates  are 
much  smoother  than  the  wildly  oscillating  estimates  of  the  pure  time 
delay  tracker,  which  are  caused  by  the  switching  of  the  semi -Markov 
weights . 

The  superior  performance  of  the  hybrid  filter  is  again  evident 
in  Figure  6.3.4  where  its  estimates  of  the  Source  vertical  position  are 
always  within  30  feet  of  the  actual  value  (400  feet) .  Since  this  error 
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Produced  by  Both  Fillers 


Estimates  of  Source  Radial  Velocity  Produced  by  Both  Filters 


Figure  4.3.4  Estimates  of  Source  Vertical  Position  Produced  by  Both  Filters 


of  30  feet  is  comparable  to  the  physical  height  of  a  submarine,  the 
error  is  entirely  acceptable.  Contrasting  with  this  small  error  are 
the  relatively  large  errors  produced  by  the  pure  time  delay  filter, 
which  for  values  of  t  _>  15  minutes  are  in  excess  of  100  feet. 

Note  how  in  Figures  6.3.2  and  6.3.4  (and  to  a  lesser  extent  in 
Figure  6.3.3)  the  estimates  provided  by  the  pure  time  delay  filter  are 
progressively  getting  poorer  and  poorer.  For  example,  as  time  increases, 
the  amplitude  with  which  these  estimates  oscillate  is  getting  steadily 
larger;  near  t  *  0,  the  estimates  are  relatively  close  to  the  true  values 
but  with  each  iteration  the  estimates  diverge  more  and  more.  The  under¬ 
lying  cause  of  this  divergence  is  that  the  signal- to-noise  ratio  on 
is  rapidly  deteriorating  to  the  point  where  the  noisy  measurements, 

z  ,  are  virtually  useless.  For  example,  over  the  period  t  *  0  to  t  * 

1 

25  minutes,  decreases  from  an  initial  value  of  18  m.sec  to  a  final 
value  of  6  m.sec.  With  additive  white  Gaussian  measurement  noise  having 
a  5  m.sec  standard  deviation,  an  initially  marginal  quality  for 
rapidly  deteriorates  to  the  point  where  the  additive  measurement  noise 
is  of  the  order  of  the  quantity  being  measured.  The  following  specific 
examples  taken  from  the  computer  simulations  used  to  produce  Figures 
6.3.2  -  6.3.4  serve  to  illustrate  this  point. 

t  (Actual  Value  in  m.sec)  z  (Measured  Value  in  m.sec) 

T1 

16.3  23.3 

13.7  19.6 

10.7  16.6 

7.0  0.1 

The  hybrid  filter,  on  the  other  hand,  performs  remarkably  well  when  con¬ 
fronted  with  these  same  measurements.  The  progressive  deterioration  in 
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the  quality  of  the  estimates  encountered  with  the  pure  time  delay  filter 
is  noticeably  absent  from  the  hybrid  filter's  estimates.  In  fact,  the 
high  quality  of  the  estimates  produced  by  the  hybrid  filter  at  the  be¬ 
ginning  of  these  plots  when  the  z  measurements  are  of  marginal  quality, 

T1 

is  maintained  throughout  the  remainder  of  the  plots  even  when  z  becomes 

T1 

virtually  useless. 

The  reason  for  the  superior  performance  of  the  hybrid  filter  is 
obvious;  the  high  quality  of  the  estimates  (6.2.1)  produced  by  the  Adap¬ 
tive  Doppler  tracking  filter  (which  i3  immune  to  the  debilitating  effects 
of  poor  quality  SNR  on  the  time  delay  measurement  data)  is  sufficient  to 
arrest  the  deterioration  caused  by  the  decreasing  SNR  on  t^.  Thus  the 
initial  high  quality  of  the  estimates  produced  by  the  reduced  order 
sonar  time  delay  filter  is  maintained  throughout  the  entire  time  period. 

A  factor  which  aids  both  filters  in  the  above  tests  is  the 
high  SNR  of  Che  measurements.  The  ocean  depth  of  6000'  produces 
values  of  in  the  range  721  m.sec  -  286  m.sec. 

In  the  next  series  of  tests  the  ocean  depth  is  reduced  to  3000' 
with  all  initial  conditions  and  filter  parameters  remaining  unchanged. 


This  reduction  in  the  depth  of  the  ocean  reduces  the  values  of 
to  one-fifth  of  their  previous  values.  The  values  of  remain 
the  same  as  before.  Thus,  coupled  with  an  already  very  low  SNR  on 
is  a  realtively  low  SNR  on  t^.  In  addition,  the  Source  executes 
e  maneuver  which  causes  its  radial  velocity  (with  respect  to  the  moving 
Observer)  to  change  from  an  initial  16  ft/sec  to  a  final  10  ft/sec. 
These  conditions  collectively  provide  a  much  more  rigorous  test  of 


both  filters 
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Figure  6.3.7  Estimates  of  the  Source  Radial  Velocity  Produced  by  Both  Filters 

with  Very  Low  SNR  on  i,  and  t 


Figure  6.3.9  Estimates  of  Source  Vertical  Position  Produced  by  Both  Filters 

with  Very  Low  SNR  on  t,  and  x_ 


6.3.10  Continuation  of  Figure 


The  effect  of  the  very  low  SNR  on  t2  in  addition  to  that  on 
is  immediately  evident  from  Figure  6.3.5  and  its  continuation  Figure 
6.3.6  where  the  percent  error  produced  by  the  pure  time  delay  filter  is 
from  three  to  five  times  greater  than  the  corresponding  values  in  Fig¬ 
ure  6.3.2.  These  large  percent  errors  (132  -  252)  stand  in  marked  con¬ 
trast  to  the  very  small  percent  error  (02  -  32)  produced  by  the  hybrid 
filter  when  presented  with  the  same  measurement  data. 

In  Figures  6.3.7  and  6.3.8  the  hybrid  filter's  estimates  of  the 
Source  radial  velocity  which  are  generally  within  ±1  ft/sec  of  the  true 
value  are  far  superior  to  the  highly  oscillatory  estimates  produced  by 
the  pure  time  delay  filter.  The  hybrid  filter  also  responds  much  more 
rapidly  to  the  Source  maneuver  at  t  *  4  minutes.  For  example,  by  t  *  6 
minutes,  its  estimates  have  converged  to  che  new  value  of  the  Source 
radial  velocity  whereas  the  other  filter  exhibits  a  very  sluggish  re¬ 
sponse  which  lasts  until  t  *  12.5  minutes.  Thus  the  response  time  of 
the  hybrid  filter  is  one- fourth  that  of  the  pure  time  delay  filter. 

This  fast  response  by  the  hybrid  filter  is  a  direct  result  of  embedding 
its  adaptive  feature  in  the  Doppler  frequency  filter. 

The  estimates  of  the  Source  vertical  position  produced  by  the 
hybrid  filter  in  Figures  6.3.9  and  6.3.10  continue  to  be  generally 
within  30  feet  of  the  true  value.  This  error  is  entirely  acceptable  for 
the  reasons  stated  earlier.  The  estimates  produced  by  the  pure  time 
delay  filter,  on  the  other  hand,  are  in  error  by  several  hundred  feet 
for  values  of  t  ^  17  minutes. 

The  progressive  deterioration  of  the  pure  time  delay  estimates, 

originally  encountered  in  Figures  6.3.2  -  6.3.4,  is  even  more  pronounced 

in  Figures  6.3.5  -  6.3.10.  The  reason  for  this,  of  course,  is  that  in 
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addition  to  a  very  low  SNR  on  the  SNR  on  t ^  is  much  lower  in  the 
latter  set  of  figures  than  is  the  case  in  Figures  6.3.2  -  6.3.4.  The 
immunity  to  these  low  SNR  ratios  which  the  Adaptive  Doppler  frequency 
filter  enjoys  is  what  enables  the  hybrid  filter  to  maintain  its  consis¬ 
tently  superior  performance  in  both  series  of  tests. 

These  high  quality  estimates  of  both  the  Source  radial  and  ver¬ 
tical  positions  produced  by  the  hybrid  filter  under  very  low  SNR  on 
and  ?2  tend  to  obscure  the  fact  that  the  reduced  order  sonar  time  delay 
filter,  like  the  pure  time  delay  filter,  is  also  an  Extended  Kalman  fil¬ 
ter.  It  is  well  known  that  the  Extended  Kalman  filter  algorithm  dis¬ 
plays  large  errors  and  biases  under  low  signal-to-noise  ratios  . 

This  fact  makes  the  performance  of  the  hybrid  filter  in  this  chapter 
all  the  more  remarkable.  In  Figures  6.3.5,  6.3.6,  6.3.9  and  6.3.10, 
with  a  low  SNR  on  Tj  and  an  SNR  on  which  makes  z ^  virtually  use¬ 
less  most  of  the  time,  the  estimates  produced  by  the  reduced  order  sonar 
time  delay  filter  show  a  relatively  small  degradation  compared  to  that 
suffered  by  the  pure  time  delay  filter.  In  addition,  this  good  per¬ 
formance  by  the  reduced  order  filter  is  maintained  in  the  presence  of 
an  abrupt  Source  maneuver. 

6.4  Conclusion 

In  this  chapter  an  integrated  Adaptive  Doppler  frequency- 
reduced  order  sonar  time  delay  tracking  filter  is  developed  by  com¬ 
bining  the  filters  developed  in  Chapters  3  and  5.  This  hybrid  filter 
exploits  the  best  features  of  both  constituent  filters  by  using  the 
high  quality  estimates  of  the  Source  radial  velocity  and  control  input 
provided  by  the  Adaptive  Doppler  frequency  filter  to  improve  the 
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estimates  of  the  Source  radial  and  vertical  position  produced  by  the 


pure  time  delay  filter.  For  example  the  immunity  to  sonar  time  delay 
measurement  noise  which  the  Adaptive  Doppler  frequency  filter  enjoys, 
has  carried  over  to  the  reduced  order  sonar  time  delay  filter.  Under 
low  SNR  conditions  on  both  time  delays,  the  large  errors  in  the  es¬ 
timates  of  the  Source  radial  and  vertical  positions  produced  by  the 
pure  time  delay  filter  of  Chapter  3  are  noticeably  absent  with  the 
reduced  order  filter.  In  addition,  the  hybrid  filter's  response  time 
to  a  maneuver  by  the  Source  is  found  to  be  one-fourth  that  of  the  pure 
time  delay  filter. 

The  parameter  values  used  in  this  study  of  the  hybrid  filter 
may  not  be  such  as  to  yield  optimal  filter  performance;  therefore,  any 
further  investigation  of  this  filter's  characteristics  should  include 
a  sensitivity  analysis  of  these  various  parameters.  Such  an  analysis 
might  well  yield  a  more  "optimal"  set  than  the  one  used  here. 


SUMMARY 

Modeling  of  the  target  control  input  as  a  series  of  partially 
overlapping  Gaussian  curves  worked  well  against  air  targets  in  the  past, 
it  was  decided  to  apply  this  technique  to  the  problem  of  passively 
tracking  an  underwater  maneuvering  Source  from  a  moving  Observer.  In 
order  to  do  this,  it  was  necessary  to  carefully  evaluate  the  benefits 
that  the  various  coordinate  systems  have  to  offer  given  the  generally 
higher  nonlinear  types  of  measurement  data  available  for  passive 
tracking.  When  polar  coordinates  were  chosen  as  the  most  suitable 
coordinate  system,  a  linearized  polar  model  for  the  maneuvering  Observer- 
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Source  scenario  was  developed  which  solved  the  difficult  problem  of 
compensating  for  the  Observer's  own  motion.  This  choice  of  polar  co¬ 
ordinates  represents  a  significant  departure  from  the  existing  trends 
in  the  literature  which  generally  show  an  unquestioned  preference  for 
rectangular  coordinates.  Indeed,  one  of  the  major  contributions  of 
this  report  is  to  show  that  for  the  small  added  complexity  needed  to 
compensate  for  Observer  motion  in  polar  coordinates,  considerable  fil¬ 
ter  simplification  can  be  bought  in  the  processing  of  all  the  types  of 
available  data — sonar  time  delays,  bearing,  and  Source-emitted  fre¬ 
quency  spectra.  With  rectangular  coordinates,  on  the  other  hand,  where 
compensation  for  Observer  motion  is  extremely  simple,  the  filter  is 
generally  more  complex  than  in  the  polar  case.  It  is  shown,  for  ex¬ 
ample,  that  in  processing  sonar  time  delays  the  linearized  measurement 
matrix  using  rectangular  coordinates  is  at  least  one  order  higher  than 
that  needed  with  polar  coordinates. 

Using  partially  overlapping  Gaussian  curves  to  model  the  Source 
control  input  and  with  measurement  data  consisting  of  both  surface  and 
bottom  reflected  sonar  time  delays,  the  estimates  of  the  Source  radial 
velocity  produced  by  the  polar  filter  were  found  to  exhibit  some  os¬ 
cillations  even  under  high  measurement  SNR  conditions. 

As  a  result  of  these  oscillations,  a  new  technique  using  the 
Doppler  effect  was  developed  whereby  high  quality  radial  velocity  esti¬ 
mates  are  obtained  from  the  frequency  spectra  emitted  by  the  Source. 

The  processing  of  these  spectra,  which  are  either  unknown  constant  ton&i, 

OX  naxxowband  Xandom  pX0ce.4A£4,  is  a  highly  nonlinear  estimation  problem, 
which  when  implemented  on  the  polar  filter  became  linear  and  time  varying. 
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Another  important  result  of  this  report  was  to  show  how  high  quality 
information  about  the  Source  radial  velocity  and  control  input  can  be 
obtained  from  these  unknown  spectra  using  the  Doppler  effect.  In 
addition,  and  in  spite  of  the  nonlinearity  of  these  noisy  frequency 
measurements,  the  information  was  obtained  using  only  linear  estimation 
techniques . 

Next,  an  integrated  filter  was  developed  by  combining  the 
Doppler  frequency  filter  with  a  reduced  order  sonar  time  delay  filter. 
The  estimates  produced  by  this  integrated  filter  are  generally  superior 
to  those  produced  by  the  pure  time  delay  filter.  This  superior  per¬ 
formance  was  maintained  under  low  SNR  on  both  sonar  time  delays  and 
represents  a  major  contribution. 
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